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I. INTEODOCTION 



An overall rise in fuel prices has led to an increasing 
interest in the design of autopilots for ships. The purpose 
of the automatic steering control is to minimize the propul- 
sion losses, vhich are caused by added drag due to steering 
of the ship. iMinimizing a performance criterion based on 
added drag due to steering can reduce fuel consumption. 
Claims by many researchers indicate that a carefully 
designed controller could save from one to two percent of 
fuel. For large containerships this could amount to more 
than ^100,000.00 per jear savings. 

To study the optimization problem, models of both the 
ship and its operating environment are required. What type 
of computer model should be used to represent the ship? 
Chapter two addresses the develoiJmeat of several models. 
Since the best model was desired it was decided to use the 
equations of motion to simulate the ship in our Fortran 
program. The basic Nomoto models give an adequate descrip- 
tion of ship steering dynamics for design. The Ncmoto 
second- and third-order models were developed from the equa- 
tions of motion as defined by a series expansion including 
all terms (both linear and nonlinear) for which hydrodynamic 
coefficients were available. An interactive program that 
utilize'! the Nomoto mcdels to model the ship was also used. 
Two independent programs were developed to aid in the design 
of the controller. 

ihat is an adequate cost function which represents the 
added drag due to steering? Chapter three addresses the 
classical cost function used by many researchers. 

Since a variety cf control algorithms are possible one 
must ask if one algorithm provides a lower minimal cost than 
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another. Chapters fear, five and six address the selection 
of the controller which provides the minimum value of added 
drag due to steering. 

Ship dynamics change with operating conditions such as 
ship speed , sea state , and encounter angle. Therefore an 
adaptive controller must be used to provide minimum added 
drag due to steering. Chapter seven development of an 
approach to an adaptive controller utilizing satellite 
inf or ma ticn- 

Conclusions were drawn from these experiments and are 
presented in Chapter eight. This thesis investigated only 
course keeping with emphasis on minimizing rudder and yawing 
activity to reduce fuel consumption. Presented in this 
Chapter are recommendations for future study where the 
objective is track fcllowing which would be important for 
ships required to follow stringent routes. It is also impor- 
tant for other systems such as satellites, missiles, 
aircraft, where the controller minimizes yaw error to keep 
the system on track. 
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II, COMPOTER MODELS 



The Bodel vhich test represents ship-steering dynamics 

is a Taylor’s series expansion of the force and moment rela- 

% 

tionships around a selected steady-state operating jJoint. 
The resulting eguaticns are commonly known as the equations 
of motion [Bef* 1], A computer program was developed using 
known available data on the hydrodynamic coefficients for 
the Sl-7 containership to provide a computer simulation of 
the ship. The computer program is shown in Appendix A. 
figure 2.1 shows the block diagram. Small yaw command 
angles are used, for example Yawc= 1.0 / 57.296 represents a 
yaw command change of one degree. 



I 




Figure 2.1 BLOCK DIA6BAH 
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To obtain the Noaoto second- and third-order transfer 
functions from the equations of motion, the function mini- 
mization subroutine vas used to obtain the coefficients. 
Figure 2.2 shows the scheme used to obtain the Ncmoto 
transfer functions. The computer program is shown in 
Ippendix A. 




The Boaoto models were checked against analytic results 
from linearized equations. 

Proceeding to the second-order Homoto equation: 
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yl^[S)/d (S) = K/S* (1 +I*S) 



( 2 . 1 ) 



Deriving the second-ctder Nomoto transfer function froic the 
ya¥ equation only, the result is 

^(S)/5(S) = 0. 040893/S* ( 1 + 8.539932*3) 
and using function minimization as in Figure 2.2 
yp[S)/d{S) = 0.0409221/S* (1 + 8.5520782*5) 
and the agreement is obvious. Osing function minimization 
with both yaw and sway equations with linear terms only, the 
results are: 

^(S)/a(S) = 0. 1072741/S* (1 + 3 1.9 199524*S) 

If the nonlinear terms are included but the perturbation is 
small 

^(S)/a(S) = 0. 1072082/S* (1+31.8907013*3) 
and it is clear that the nonlinear terms contribute little. 
Proceeding to the third-order Nomoto equation: 

tA(S)/5(S) = K* {1 + IZ*S)/S* (1+TP1*S) * (1+TP2*S) (2.2) 

Ihe parameters were calculated and checked by using function 
minimization as in Figure 2.2. The results are given in 
Table 1. It is clear that the answers obtained by function 
minimization agree closely with the analytic solutions. 



TABLE 1 

THIED-OEDEB NOMOTO MODEL FOE THE SL-7 



speed K TZ TP1 TP2 

knots calc comp calc comp calc comp calc comp 

16 .0738 .0738 22.57 22.95 12.946 12.946 107.583 107.583 

23 .1067 .1061 15.67 15.70 9.014 9.006 75.130 74.846 

32 .1477 .1477 11.28 11.28 6.470 6.467 53.793 53,793 
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Analytical equations used to calculate second-order 
Nomotc transfer function coefficients are: 

K =N /N ; 1 =- (I -N . ) /N 

6 r z r ' r 

Analytical equations used to calculate third-order 
transfer function coefficients are: 

K =(N^-U )/(N -N * {Y -H + U) /Y ) 

6 v6 v rv r v 

TZ =-( (E-Y*) -N-*Y. ) /(Y ♦N. -H *Y. ) 

TP1*IP2=-((M-Y*)*{1^”N*)-N‘*Y-)/(N *(Y -M*U)-Y ♦N ) 

' z r V r v r v r 

TP 1 + 1P2= ( {M-Y •) + (I -N-)*Y +N. *(Y -a*U)+Y**N ) 

v'r zr vv'r rv 

/(N «(Y -M*U)-Y *N ) 

The nondimensionalized hydrodynamic coefficients for the 
SL-7 ccntainership are shown in Table 2. 



TABLE 2 

HYDBODYNAHIC COEFFICIENTS FOB THE SL-7 



axial 


force 


lateral 


force 


moment 


2 -axis 


A ^ 




-0.0001 


Y* 

V 


-0.00758 


N* 

V 


-0.00213 


1 » 

uu 


SI 


-0,0003 


Y* 

r 


0.0023 


N» 

r 


-0.00105 


X* 

vr 


= 


0.0039 


Y* = 

6 


0.00145 


"'6 ' 


-0.0007 


X • 

^ vv 




-0.0012 


Y* 

vvr 


0.01 


N’ 

vvr 


-0. 015 


X ’ 

^ 66 




-0.0005 


Y» 

vrr 


-0.008 


B* 

vrr 


-0. 008 








Y f 

vvv 


-0. 03 


N* = 

vvv 


0, 0 1 


\ 






Y t = 

rrr 


0.003 


N* 

rrr 


-0,006 








^ 666 


-0.0005 


N * = 

66 6 


0.0001 



16 



III. COST FDN CTIQN 



In recent years, many have studied the problem of 
[fief. 2] [fief. 3] [fief. 4] [fief. 5] [fief. 6 ] [fief. 7] 
[fief. 8 ] [Ref. 9] [fief. 10] [Ref- 11] [fief- 12] [fief- 13] 
[fief. 14] optimizing an automatic ship-steering controller 
for minimum fuel consumption. It is well known that addi- 
tional drag is introduced by steering and that both the 
rudder motion and the yawing motion contribute to this added 
drag- A measure of the added drag given as a cost function 
is 



T 



VT/( 

' O 


X ♦ 2 4- 52 ) (it 


(3.1) 


where 


= yaw error 






6 - rudder angle 






X = weighting factor 





While this expression is an approximation, it is conven- 
ient for shipboard use because and 5 are readily measur- 
able. There is no general agreement on numerical values for 
the weighting factor, X , and in this study the values used 
were chosed from the work of fi.E. Reid [Ref. 7] for the 
SL-7- 

The weighting factors for the operating range of the 
ship are shown in Table 3- 
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TABLE 3 

HEIGHTIHG FACTOfi 



ship speed 



(knot s) 



weighting factor 



16 

23 

32 



16.796 

8.128 

4.2 



Reid’s work shows the relationship of weighting factor 
to the closed-loop natural freguency, mass, pivot point, 
ship speed, ^*66 i^ydrodynaraic coefficients. It is 

shown in Equation 3.2. Reid chose a closed- loop natural 
freguency of 0-05 rad/sec which experimentally shewed at 
this frequency, the ^«eighting factor in the cost function, 
provided good representation of the added drag due to 
steering. 

X = 2*f5*(1 + X” )*(CP/L)*u>2 /(p/2) * (1-^X,, *02) (3.2) 

VIT 0 0 
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IV. COHTBOLLEH DE SIGN OSIHG ICSOS 



The Interactive Control System Optimization and 
Simulation (ICSOS) package finds optimum values for unknown 
(free) parameters in a control system design problem and/or 
performs simulation of the system- An example of usage of 
ICSOS is shown in Appendix B. 

In preliminary studies ICSOS was used with Nomoto models 
to study controller characteristics in calm water. The func- 
tion minimization subroutine adjusted the controller parame- 
ters to ninimize the cost function. Figure 4-1 shows the 
scheme used to evaluate the controller. parameters. 




Pigure 4.1 OPTIEIZATIOB OF CONTBOLLEB 



/ 
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Reid [Ref* 7] uses the second-order Nomoto model of 
equation 2.1 for the SL-7 and also uses a ccntrcller 
descriled by 

Gc (£) = K1*(1+T1*S) / (UT2*S) (4.1) 

His results are given in Table 4. 



TABLE 4 

REID'S BESOLTS 



speed 


plant 


weighting 


controller 


gai ns 


knots 


K T 


factor 


K1 


T1 


T2 


16 


0. 1084 90.36 


16.796 


0.4556 


89.51 


10.06 


23 


0.1556 64.67 


8.128 


0.3769 


62.60 


8.308 


32 


0.2167 45.45 


4-2 


0-3188 


44.92 


7.066 


Using 


this plant and 


weighting 


factor 


values 


but apply 



ICSOS, results were obtained and shown on Table 5. 



TABLE 5 
ICSOS BESOLTS 



speed 


plant 


weighting controller 


gains 


cos t 


knots 


K T 


factor K1 


T1 


T2 


J min 


16 


. 1084 90.36 


16.796 .454616 


90.3459 


10.0215 


340- 864 


23 


.1556 64.67 


8.128 .373171 


64-6658 


8-4640 


139.9916 


32 


.2167 45.45 


4.2 .318645 


45.4475 


7.0662 


60.828 



In each case the controller zero (1/T1) cancels the 
plant pole (1/T). Additional experiments consisting of 
inserting arbitrary numbers in the Nomoto equation and 
repeating the computer run indicated that this will always 
be true. That is, to minimize the cost the plant pole is 
cancelled and a new pole location determined with appropri- 
ately adjusted gain. 
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The sittple ccntrcller of Eguation 4.1 is an arbitrarily 
chosen structure. To determine the effects of more complex 
controllers three additional structures were chosen as shown 
in Figure 4.2. Each of these was used with the Ncmoto 
second-order model for the ship at each of the indicated 
speeds. The results are shown in Tables 6, 1, and 8. 




Figure 4.2 VIBIOOS STBOCTOBES FOB COBTBOLLEBS 

These results are very interesting, at 16 knots the 
controller gain (K1) , controller zero (1/T1) and controller 
pole (1/T2) cire essentially the same for all structures. For 
structure E, which includes an additional pole, the function 
aiinimization subroutine tries to drive the additional pole 
to infinity, and no doubt would have done so if the calcula- 
tions had continued. For structure C, which has two poles 
and two zeros, a zero and pole cancel indicating that they 
are not needed or wanted. For structure D, the integrator 
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TABLE 6 

SIflULATION EE5ULTS - STEADY STATE 600 SECS 



CALH HATEB FOB VABIOOS COHTBOLLEBS 
FOE FIXED SHIP SPEED { 16 KNOTS ) 
HOHOTO SECONE-OEDEB HODEL (K=. 1084 , T=90. 36) 
X= 16.796 , OPTIMAL PAEAMETEB GAINS OF 
VABIOOS CCNTBOLLEES , COST FUNCTION 



contr 


controller 


gains 




cost. 


K1 


T1 


T2 


T3 T4 


Ti 


J min 


A .454616 


90.4359 


10.0215 


— — 


— 


340-864 


B -444101 


90.2950 


9.8566 


. 0 1 


— 


341.046 


C .454511 


90-3685 


10.0224 


23.085 23.084 


— 


340.864 


D .454581 


90.3719 


10.0222 


— — 


1E09 


340-864 



TABLE 7 

SIflOLATION BESOLTS - STEADY STATE 600 SECS 



CALM HATifi FOB VABIOOS CONTBOLLEBS 
FOB FIXED SHIP SPEED { 23 KNOTS ) 
NOHOTO SECONE-OBDEB MODEL (K=. 1556 . T=64. 67) 
X= 8.128 , OPTIMAL PABAMETEB GAINS OF 
VABIOOS CCNTBOLLEES , COST FUNCTION 



con tr 

K1 

A .373171 
B .3«l0024 
C .373139 



ccntrcller gains 
T1 T2 



64.66579 

79.65872 

64.66855 



8-463957 
8.889204 
8. 463497 



T3 T4 

.01 

25.9719 25.9738 



cost 
J min 

139.9916 

140.9338 

139.9991 



TABLE 8 

SIMULATION BESOLTS - STEADY STATE 600 SECS 



CALH BATIB FOB VABIOOS COHTBOLLEBS 
FOB FIXED SHIP SPEED ( 32 KNOTS ) 
NOflOTO SECOND-OBDEB MODEL (K=. 2167-T=45. 45) 
X= 4.2 . OPTIMAL PABAMETEB GAINS OF 

VABIOOS CCNTBOLLEES , COST FONCTION 



contr 

K1 


controller gains 
TI T2 T3 


T4 


cost 
J min 


A .318645 


45- 44747 


7-06617 




60.828 


B .318 


45.45 


7.066 .05 


— 


60.933 


C .318678 


45.57511 


7.06790 50-1829 


50.04832 


60.828 



f 
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gain is driven to zero. The same pattern of results is 
obtained at 23 knots and 32 knots. Note that in all cases 
the irinimum cost is essentially the same, as vould be 
expected since all controllers are the same. 

Using the computer method of Figure 4.1 and the Ncmoto 
third-order models of Table 1, controllers A, B, C of Figure 
4.2 were optimized.. The results are shown in Tables 9, 10, 
and 11. 



TABLE 9 

SIflOLATION RESULTS - STEADY STATE 600 SECS 



CALM lATEH FOR YARIOOS CONTROLLERS 
FOE FIXED SHIP SPEED ( 16 KNOTS ) 

NOMCTO THIRD-ORDER MODEL 
(K=.073812,TZ=22.5673,TP1=12.9458,TP2=107.5853) 
16.796 , OPTIMAL PARAMETER GAINS OF 
VARIOUS CCNTROLLEES , COST FUNCTION 



contr 


controller gains 








K1 


T 1 


T2 


T3 


T4 


A 


0.6446104 


90.CS94 


15-27712 






B 


0.6441367 


84- 826 


15.78691 


-24598 


- 


C 


0.6151139 


107. 5782 


8. 73520 


12.9368 


24.9676 



cost 
J min 

370.4023 
374.3808 
369- 9297 



TABLE 10 

SIMULATION RESULTS - STEADY STATE 600 SECS 



CALM HATER FOR VARIOUS CONTROLLERS 
FOB FIXED SHIP SPEED ( 23 KNOTS ) 

NOMCTO THIRD-ORDER MODEL 
(K=, 1067.TZ=15. 675,TP1=9.014, TP2=75. 13) 

X= 8.128 , OPTIMAL PARAMETER GAINS OF 
VARIOUS CCNTROLLEES , COST FUNCTION 

contr ccntroller gains cost 

K1 T1 T2 T3 T4 J min 

A 0.5224258 63.13609 12.72212 - - 152.2920 

B 0.5216467 64.93709 12.63218 .0505174 - 152.5333 

C C. 5001907 75-14852 6.527490 9.039928 18.260 152.2800 



Of iiajcr interest is the fact that the difference in 
"cost” between A, B, C is less than one per cent. At each 
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TABLE 11 

SIMULATION HISOLIS - STEADY STATE 600 SECS 



CALM HATEB POE 7ABIOOS CONTfiOLLEBS 
POE PIXEC SHIP SPEED { 32 KNOTS ) 

NOMCTO THIED-OBDEB MODEL 
(K=. 1477 1 -TZ = 11. 2833, TP1=6. 4699, TP 2=53. 7931) 
A. = 4.2 , OPTIMAL PABAMETEfi GAINS OP 

VARIOUS CCNTEOLLEBS , COST FUNCTION 



contr 



A 

B 

C 



contrcller gains 
K1 T1 T2 

0.427633 48.66C48 10. 74485 
0.298732 89.40696 15-01033 
0.417991 53-69961 4,970016 



cost 

T3 T4 J min 

68.09039 

.0597786 - 69.32355 

6.294354 13.85724 68.04735 



speed (16,23,32 knots) controller C is "BEST”, tut the 
difference is slight. Examining the parameter values 
obtained fcr contrcller C, it is seen that at all three 
speeds both poles of the ship are essentially cancelled by 
zeros of the controller. 

These results seem to indicate that the dynamics cf the 
plant determines the optimum structure for the contrcller. 

Using a state-feedback controller and Nomoto third-order 
models of Table 1, the contrcller was optimized for various 
ship speeds. Figure 4.3 shows the scheme used to evaluate 
the state-f eedback controller. 

Using the scheme of Figure 2.2, with no change in cost 
function or weighting, the optimal gains and costs were 
determined as shown in Table 12, 

Hhen comparing the state-feedback controller with 
contrcller C, it is seen that at each speed contrcller C has 
a lower cost. Among the controllers consired, controller C 
is "BEST" when using the Nomoto third-order model, although 
the differences in cost are not dramatic. 
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Figure 4.3 OPTIBIZATIOH OF STATE FEEDBACK COHTfiOLLEfi 



TABLE 12 

SIflOLATION BESOLTS - STEADY STATE 600 SECS 



CALfl VATEfi FOB VABIOOS SHIP SPEEDS, 

OPTIBAL PABABETEB GAIHS FOB 
STATE-FEEDBACK COHTBOLLEB 

speed Hcmoto third-crder plant weighting controller cost 
knots K TZ TPl TP2 factor K1 K2 J ain 

16 .0738 22.567 12-946 107.583 16.796 4.426 78.004 370.711 

23 .1067 15.675 9.014 75.13 8.128 3.103 45.649 152.596 

32 .1477 11.283 6.470 53.793 4.2 2.240 27.896 68-2513 
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V. COMTROLLI B DESIG H DSIHG FORTRA H PROGRAM 



The Fortran progxaa referenced in Chapter two which 
provided a computer simulation of the SL-7 ship was modi- 
fied* A function minimization subroutine was coulped to the 
simulation and used the subroutine to adjust controller 
parameters to minimize the cost function and to evaluate the 
minimum cost* Figure 5* 1 shows the scheme used to evaluate 
the controller parameters* This program was used for compar- 
ison to ICSOS* The computer program is shown in Appendix A* 

I 1 




Figure 5.1 OPTIMIZATION OF CONTROLLER OSING FORTRAN PROGRAM 
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Using the computer method of Figure 5. 1 and the nonli- 
near equations of motion, controllers k, B, C of Figure 4.2 
were optimized. The results are shown in Tables 13, 14, and 
15. 



TABLE 13 

SIflOLATION RESULTS - STEADY STATE 600 SECS 



CALM HATER FOR VARIOUS CONTROLLERS 
FOR FIXED SHIP SPEED ( 16 KNOTS ) 
ECOATIONS OF MOTION 
16.796 , OPTIMAL PARAMETER GAINS OF 
VARIOUS CCNTROLLERS , COST FUNCTION 

contr controller gains 

K1 T1 T2 T3 T4 

A .646401 89.81704 15.381699 

B .620050 90-67294 15.542297 0.9201336 

C -617326 107.1494 8.597198 13-353928 25.21362 



TABLE 14 

SIMULATION RESULTS - STEADY STATE 600 SECS 



CALM HATER FOR VARIOUS CONTROLLERS 
FOR FIXED SHIP SPEED ( 23 KNOTS ) 
ECOATIONS OF MOTION 

X- 8.128 , OPTIMAL PARAMETER GAINS OF 
VARIOUS CCNTROLLERS , COST FUNCTION 

contr controller gains 

K1 T1 T2 T3 T4 

A .522106 66.33122 12-83327 

B .4S5869 66.15152 13.01183 0-92783 

C .503967 74.79771 6.65880 9.20533 18-4022064 



These results agree with those obtained by ICSOS and 
controller C provides the minimum cost. 

If the assumption that the steering dynamics of the ship 
is adequately modeled as a second-order system is valid, 
then only two states are needed for feedback. For a third- 
order system three states are required. The controller 
structures are shown on Figure 5.2 and 5-3. Using the scheme 



cost 
J min 

0-4640879 

0-4857854 

0-4636095 



cost 
J min 

1. 128189 
1. 173323 
1. 126307 
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TABLE 15 

SIHOLATIOM EESOLTS - STEADY STATE 600 SECS 



CALB BATEB FOR YABIOOS COHTBOLLEBS 
FOB FIXED SHIP SPEED { 32 KHOTS ) 

EQOATIOHS OF HOTIOH 

X* 4.2 . OPTIMAL PABABETEB GAIMS OF 

FAHIOOS CCHTfiOLLEBS , COST FONCTIOH 

contr controller gains cost 

K1 T1 T2 T3 T4 J min^ 

A .428404 48.65540 10.814426 - - 0-2072417 

B -298732 89.40696 15-010330 0-01 - 0.2118334 

C .417333 53.09654 5.096548 6-474857 14.0205 0.2071124 



of Figure 5.1 with no change in cost function or weighting, 
the optimal gains and cost were determined as shown in 
Tables 16 and 17. Comparing costs, there is little differ- 
ence between the twc state system and the three state 
system. The comparison between state feedback with 
controller C, it is seen that at each speed controller C is 
better, but not much better. 




Figure 5.2 TVO STATE SXSTEfl 
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Pigure 5*3 THEBE STATE SISTEfl 



TABLE 16 

; 

SIHOLATIOH EISULTS - STEADY STATE 600 SECS 

CALH VATEB FOB TABIOOS SHIP SPEEDS 
EQUATIONS OP MOTION 
OPTIBAL PA BA METE B GAINS FOB 
TIO STATE SYSTEM 



speed 

knots 


gains 

K1 K2 


weighting 

factor 


cost 
J min 


16 


4.4033689 


77.5041656 


16.796 


1.128771 


23 


3.0889006 


45.2637787 


8.128 


. 4646050 


32 


2.2342062 


27.6808014 


4.2 


.2075207 



/ 



Note that for the Nomoto model studies yaw error and 
rudder angles were measured in degrees; when the equations 
of notion were simulated yaw error and rudder were in 
^ radians. Thus the numerical values of the cost, J, are 
different. y 

Transient response plots for controllers A, B, C, and 
three state-feedback at ship speed 32 knots are shown in 
Figures 5. 4 - 5.9. 
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TABLE 17 

SIflOLATION BESOITS - STEADY STATE 600 SECS 



CALfl VATEfi FOR 7ABI00S SHIP SPEEDS 
ECOATIONS OF HOTIOM 
lAL 



OPTIBAL PABASETEB GAIHS FOB 
THREE STATE SYSTEM 



speed 

knots 

16 

23 

32 



K1 

4,8617249 

3.6630983 

2.5967150 



gains 

K2 

87.7073364 

56.2784882 

33-7511444 



K3 

99.9802704 

88-5913391 

41-3186035 



weighting 

factor 

16.796 

8.128 

4.2 



cost 
J sin 

1-128289 

-4643548 

.2074225 




Figure 5.4 YAH ts. .TIflE (controller A, C and state-feedhack) 
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Figure 5*5 BODOEF 'ws. TIHE (controller A, B, C and state-feedback) 
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32 KNOTS 

YfiH-RUOO£R VS TIM£»CO«P R 




Figure 5*6 Y18 AHD BUDOEB vs. TIHE (controller A) 
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32 KNOTS 

/ 

/ 

TflN-RUOCeR VS TlNE-COMPENSflTOR B 




Figure 5.7 TAB AHD BODDEB TS. TIME (controller B) 
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32 KNOTS 



/ 

YflW-RUOOER VS TIME-COttPCNSfiTOR C 




Figure 5.8 TAl IHD BODDEB vs. TIME (controller C) 
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YBH-RUOOCR 

,*a.o -1.5 - 1.0 ' *0.5 0.0 



32 KNOTS 

tflH-RUOOER VS TIHE-FEED8BCK COMPENS 





I ~ ~ 1 

Figure 5.9 lAI AMD BDDDEB vs. TIME, (state-feedback controller) 
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VI. C^IiOLLEB IN STATE 

A sea state generator was coupled to the Fortran 
program, so that the function minimization subroutine could 
be used in the presence of the sea state. The sea state 
generator was an elaborate program obtained from DTNSEBC. 
This program generates added mass and inertia values as 
functions of encounter frequency and also calculates the 
forces and moments. The forces and moments are generated 
and stored in a look up table which was coupled to the equa- 
tions of motion. Figure 6- 1 shows the scheme used to eval- 
uate the controller parameters. The computer program is 
shown in Appendix A. 

The optimal gains obtained by the calm water study of 
Chapter five were use as the initial guess in evaluating the 
optimal controller parameters in the presence of a seaway. 
For comparison, studies of the value of the cost function 
using calm water gains in sea state were obtained; then the 
function minimization subroutine was allowed to adjust 
controller parameters in the presence of several sea states 
and encounter angles. The entire study was done at a ship 
speed of 32 knots. The added mass and inertia change with 
respect to encounter frequency as shown in Figures 6.2 and 
6.3. Figure 6.3 is ncndimensionalized by dividing the added 
inertia by the mass of the displaced water and the square of 
the length between perpendiculars. To convert back to dimen- 
sionalized units of lb-ft-sec2, multiply the graph points by 
2. 581 El 2. Since the sea state is represented by irregular 
waves, the waves impinging on the ship hull contain the 
total energy density spectrum composed of many frequencies 
and the ship responds to an average value of added mass and 
inertia. The values used for this study was obtained at 
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encounter frequency cf 0.75 rad/sec from our sea state 
generator. This frequency gave us values for added mass and 
inertia representative of an average value. The energy 
density spectra for various sea states are shown in Figure 
6.4. The added mass for sway was changed from 2-6457E06 
lb-sec2/ft for calm water to 2.3043E06 Ih-sec^/ft for a 
seaway. The added moment of inertia for yaw was changed 
from 1-42E11 Ib-ft-sec^ for calm water to 1.5096E11 
lb-ft-sec2 for a seaway. All other hydrodynamic coeffi- 
cients were kept constant at calm water values. The results 
are shown in Tables 18 - 25. In certain sea states and 
encounter angles the calm water optimal gains performed well 
as shewn by calm water cost value when compared to sea state 
cost value. In most cases the function minimization subrou- 
tine found new gains with lower cost function values in 
seaway as compared to using calm water gains. In the calm 
water evaluation, the system was perturbed with a one degree 
course change, but the course change was not included in the 
seaway tests. The difference in cost values is attributed to 
the difference in operating conditions. 

Using the Proportional, Integral and Derivative (PID) 
controller Equation 6.1 with no change in cost, function , 
the function minimization subroutine was used to adjust 
controller parameters to minimize the cost function and 
evaluate the minimum cost. The results are shown in Table 
26. Bhen comparing the PID with controller A, it is seen 
that at each encounter angle, controller A is better. These 
results agree that in a seaway controller A provides the 
minimum cost. 

Table 27 shows comparison of the minimum cost function 
for controller A, controller C, and PID. The study was done 
at ship speed of 32 knots and at sea state 4. Controller A 
provides the minimum cost. 
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Figure 6.1 OPTIBIZATIOH OF COHTHOLLEB IH SEA STATE 

The optimal gains obtained in the presence of sea state 
was dene over using a simulating time of 600 seconds. The 
sea state program is designed to provide gradual increase in 
the forces and moments during an initial time interval-. 
This is done to minimize initial condition transients in the 
ship dynamics. There min' unavoidably be some transient 
effects, however, and these could affect the value of the 
cost, J, determined during the 600 seconds of simulation. To 
determine whether such initial transients had- any signifi- 
cant contribution to the value of J, additional simulation 
runs were made with the controller parameters fired at their 
optimal values. However, evaluation of the cost, J, was 
started only after 300 seconds of simulation had elapsed. 



38 



ADDED MASS (n0**6) 



ADDED MASS COEFFICIENTS 



SWAY WET SWAY ACCELERATION 




Figure 6.2 AIDED BASS vs. ENCOUNTEB PfiEQaENCY 
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NONDIMENSIONAL COEFFICIENTS 

0.00 0.07 0.14 0.21 0.28 0.35 0.42 




ENCOUNTER FREQUENCY (RAD/SEC) 



Figure 6.3 AODID INEfiTIA vs. EHCOaHTEB FBEQOENCI 
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AMPLITUDE SPECTRA (FT**2-SEC) 

3.0 6.0 9.0 12.0 15.0 



SEASTATE 4,5,&:6 




CIRCULAR FREQUENCY (RAD/SEC) 



Figure 6.4 BHEBGI OSHSITY SPSCTBOIi 
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TABLE 18 



SIMULATION BESOLTS - STEADY STATE 600 SECS 



EIXED SHIP SPEED (32 KNOTS) IN A SEAWAY 
SHIP MODEL: EQUATIONS OF MOTION 
YAWC=0*0 

SEA STATE 1 



CONTEOLLEE A 



encounter controller gains 

anqle 


sea state 
cost 


cost with 
calm water 


degree 


K1 


T1 


T2 


J min 


J 


0 


.4264037 


48.6554395 


10.814426 


.61745S-34 


. 61745S-34 


30 


1. 1561117 


29.3693695 


1. 4592390 


.2870198 


. 5128402 


60 


1.4033298 


10.6530075 


1. 1086883 


.1342071 


.2154726 


90 


.2969198 


58. 2413940 


1. 8758221 


. 1300669 


. 1565958 


120 


. 1761794 


299.999512 


30.7967834 


.05741726 


.0727727 


150 


2.8430557 


5.2826872 


.8887696 


.0219070 


.0939 400 


180 


1.621 1386 


14. 0782928 


2. 0712433 


.0051925 


.0095694 



TABLE 19 

SIMULATION BESULTS - STEADY STATE 600 SECS 



FIXED SHIP SPEED (32 KNOTS) IN A SEAWAY 
SHIP MODEL: EQUATIONS OF MOTION 
YAWC=0.0 

SEA STATE 2 



encounter 

angle 

degree 



CONTEOLLEE A 
controller gains 



K1 



0 .42840370 
30 .27997030 
60 .95575100 
90 1.3577642 
120 1.12C8973 
150 2-S777727 
180 .61420630 



T1 

48.6554395 
249.935059 
24. 3813629 
9.4S564080 
25. 4498596 
16.2154541 
.482041200 



T2 

10.814426 
19.857742 
2. 3079853 
1. 1068363 
4.0224676 
.56274300 
6.2521963 



sea state 
cost 
J min 

.61745E-34 

.04774852 

.04104504 

.02650556 

.04928402 

7.5751530 

.000124338 



cost vith 
calm water 
J 

.61745E-34 

.0886225 

- 0535879 

.0483197 

.0717524 

28.1294403 

.0002445 
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TABLE 20 

SIHUIATION BESOLTS - STEADY STATE 600 SECS 



EIXED SHIP SPEED (32 KNOTS) IN A SEAWAY 
SHIP MODEL: EQUATIONS OF MOTION 
YAHC=0-0 

SEA STATE 4 



CONTBOLLEB A 



encounter controller gains 

angle 


sea state 
cost 


cost with 
calm water 


degree 


K1 


11 


T2 


J min 


J 


0 


.4284037 


48.65540 


10.814426 


.620598E-34 


. 620598E-34 


30 


.9815440 


5.733036 


. 6999879 


.02854677 


.0395892 


60 


.6201209 


40. 8C556 


19. 606873 


.09375697 


.1032696 


90 


1.809746 


36. C1225 


6.324703 


1.5171340 


4. 1623011 


120 


5.195190 


18. 92513 


. 6999907 


9.991730 


48.970703 


150 


1. 446776 


16.89375 


. 5265408 


16-67052 


24.822098 


180 


.1000000 


1.000000 


20. 149999 


.00739631 


-0076657 



TABLE 21 

SIMULATION BESOLTS - STEADY STATE 600 SECS 



FIXED SHIP SPEED (32 KNOTS) IN A SEAWAY 
SHIP MODEL: EQUATIONS OF MOTION 
YAWC=0-0 

SEA STATE 6 

CONTBOLLEB A 



encounter controller gains 

angle 


sea state 
cost 


degree K1 


11 


T2 


J min 


0 


.4284037 


48.6553955 


10.8144264 


.50899E-32 


30 


2.9715786 


10.4721832 


. 5342450 


1.4287940 


60 


1.7228041 


8. 4014740 


.5141125 


1.5827220 


90 


1.8584366 


37.1672655 


.5792384 


4.5505371 


120 


3. 3422489 


106.722259 


. 9260592 


22. 108002 


150 


.2854474 


157.483887 


1 19. 981018 


.81100580 


180 


.6053379 


.75733550 


6.04484460 


.07365978 



cost with 
cal£D water 
J 

.50899e-32 
4.74724010 
3-42744920 
13.2757149 
94. 5497589 
1- 50448510 
.142564400 
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TABLE 22 

SIMULATION BESOLTS - STEADI STATE 600 SECS 



FIXED SHIP SPEED (32 KNOTS) IN A SEAHAX 
SHIP MODEL: EQUATIONS OF MOTION 
IAiC=0.0 

SEA STATE 1 

CONTBOILEE C 



encounter 


con troller 


gains 




angle 










degree K1 


T1 


T2 


T3 


T4 


0 


1.61345 


16. 8755 


25.0481 


47.3405 


1.73759 


30 


.957558 


12.6178 


7-32113 


43.3531 


8- 15752 


60 


.781984 


17.6475 


9-22485 


13.9438 


16-7663 


90 


.417332 


53- 0965 


5.09655 


6.47857 


14-0205 


120 


.417332 


53.0965 


5.09655 


6.47857 


14.0205 


150 


2. 13735 


18. 8265 


17-5778 


25. 1516 


21.1481 


180 


.957558 


12. 6178 


7.32113 


4 3. 3531 


8. 15752 



sea 
cost 
J min 

. 14E-33 
.324739 
. 1597 10 
. 148588 
.077918 
.03 1496 
.006566 



TABLE 23 





SIMULATION BESOLTS - 


STEADY STATE 600 SECS 




FIXED SHIP 


SPEED (32 KNOTS) 


IN A SEAWAY 






SHIP MODEL: EQUATIONS OF MOTION 










YAHC=0.0 












SEA 


STATE 2 












CONTBOILEE C 






encounter 


controller 


gains 




sea 


angle 










cost 


degree K1 


T1 


T2 


T3 


T4 


J min 


0 


.417333 


53. 0965 


5.09655 


6-47486 


14.0205 


. 18E-33 


30 


.849594 


19.9913 


9. 34138 


20. 3578 


13.7487 


.054338 


60 


.417333 


53. 0965 


5.09655 


6.47486 


14.0205 


.044536 


90 


.781984 


1 7.6475 


9.22485 


13.9438 


16.9438 


.033518 


120 


-880395 


21.4597 


10.9255 


9-24547 


11.0667 


.055292 


150 


-899999 


15.7103 


1. 11632 


41.3275 


2.29012 


10-7636 


180 


-440916 


- 093671 


17.7305 


25.2103 


5.04178 


.000125 



calm 

cost 

J 

. 45E-33 
.357733 
.203137 
. 1 43588 
.07791 8 
.081612 
.008172 



calm 

cost 

J 

. 1 8E-33 
.061184 
- 044536 
.048467 
.056288 
23.8522 
.001635 
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TABLE 24 

SIMOLATION BESOITS - STEADY STATE 600 SECS 



FIXED SHIP SPEED (32 KNOTS) IN A SEAWAY 
SHIP MODEL: EQUATIONS OF MOTION 
YAHC=0-0 

SEA STATE 4 



CONTBOLLEE C 



encounter 


controller 


gains 




sea 


calm 


angle 










cost 


cost 


degree K1 


T1 


T2 


T3 


T4 


J min 


J 


0 


1-22t<24 


70.3578 


61.3016 


10. 5467 


61.8215 


. 62E-34 


. 1 4E-33 


30 


.690573 


20. 1488 


20. 3214 


5.10369 


19.7841 


.034033 


.071978 


60 


.782547 


12.6178 


13.7713 


21.5637 


21.5637 


.093914 


.24436 9 


90 


2. 22895 


51.7744 


54.3190 


17.3522 


6.07814 


1. 57368 


2. 98305 


120 


3.72749 


85.4697 


40.6999 


8.52234 


1.35207 


10.3530 


37.3988 


150 


.417333 


53. 0965 


5.09655 


6.47486 


14.0205 


20.3956 


20.3956 


180 


.059166 


.286208 


19.5103 


46.4767 


22.3286 


.007397 


.099413 








TABLE 25 










SIMULATION EESULTS - 


STEADY STATE 600 SECS 






FIXED SHIP 


SPEED (32 KNOTS) 


IN A SEAWAY 








SHIP MODEL: EQUATIONS OF MOTION 










YAHC=0.0 














SEA 


STATE 6 














CONTEOILEfi C 








encounter 


controller 


gains 




sea 


calm 


angle 










cost 


cost 


degree K1 


T1 


T2 


T3 


T4 


J min 


J 


0 


2.33178 


52. 0881 


95.7892 


31. 6564 


11.6959 


.49E-32 


. 19E-31 


30 


2.08709 


73. 1270 


76.7193 


12. 1726 


16.4711 


2. 00375 


3.63379 


60 


2.00128 


71.6612 


77.6750 


13. 1 170 


17.0912 


2. 15428 


3. 41794 


90 


.957558 


12.6178 


13.7713 


72*0670 


15.4611 


5. 763 99 


8.80971 


120 


3. 10589 


81.8044 


38.2439 


91.2237 


9.14683 


24. 9099 


72.6716 


150 


1.51250 


70.3578 


61.3016 


35.5894 


6 1.82 15 


2.50971 


7.50022 


180 


1.52875 


1.87828 


43.4698 


49.9147 


11.2599 


.078894 


. 9 30 885 
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TABLE 26 

SIHOLATION EESUITS - STEADY STATE 600 SECS 



FIXED SHIP SPEED (32 KNOTS) IN A SEAWAY 
SHIP MCDEL: EQUATIONS OF MOTION 
YAWC=0.0 



SEA STATE 4 
PID CONTROLLER 



encounter 
angle K1 

30 .95263100 

60 .68631890 

90 2.5809155 

120 4.9198265 

150 1.3970823 



controller 

T1 



gains 



4.20720860 
12.5794449 
12. 4247589 
12.5986176 
15.7682953 



T2 



.69368610 
8. 2121658 
.77810380 
.67592390 
.51991 180 



sea state 
J min 

.02985619 
.09730512 
1.59 15950 
10. 708980 
17.427200 



(S) = K1 + K1*T1*S / (UT2*S)**2 



( 6 . 1 ) 



TABLE 27 

COMPARISON OF TEE MINIMOM COST 



encounter 

angle 

degree 



30 

60 

90 

120 

150 



SHIP SPEED (32 KNOTS) 
SEA STATE 4 



controller 

A 

J min 

.02854677 
.09375697 
1.5171340 
9. 99173CC 
16. 67052C 



controller 

C 

J min 

.034033 
.098914 
1. 57368 
10.3530 
20.3956 



controller 
PID 
J min 

.02985619 
.09730512 
1. 5915950 
10.708980 
17.427200 



The value obtained vas then doubled and compared with the 
result ot evaluating J over the full 600 seconds. Comparison 
of Table 28 with cost values in Tables 18, 19, 20, 21 shows 
only small differences. 

To obtain insight into the stochastic process of irreg- 
ular seas, a deterministic process was studied. The Fortran 



46 



TABLE 28 

EFFECTS DOE TO THANSIENT AND GSADOAL BUILD UP OF SEfi STATE 



INTEGRATION OF COST FUNCTION ( 300 TO 600 SECS) 

FIXED SHIP SPEED (32 KNOTS) IN A SEAWAY 

sea encounter ccntroller gains cost cost 

state angle K1 T1 T2 J min 2*J min 

1 60 1.4033298 10.650075 1. 1086683 .0641 122 . 12-82244 

2 60 .95575100 24.381363 2.3079853 .0199731 .0399462 

4 60 .62012090 40.805560 19.606873 .0515974 .1031948 

6 60 1.7228041 8.4014740 .51411250 -7906179 1.581236 



program «as modified to minimize the cost function in the 
presence of a regular sea. To allow comparison with 
previous work the encounter frequency of 0.75 rad/sec was 
used and scaled the amplitude of the regular sea to its 
prospective sea state. The entire study was done at a ship 
speed of 32 knots. The results are shown in Tables 29, 30, 
and 2 1. 

Table 29 shows that for regular seas the ccntroller 
parameters do not change significantly for different sea 
states; but as sea state increases, the cost value increases 
due to the increase in yaw moment and sway force on the 
ship. Tables 30 and 31 also show that the controller parame- 
ters do not change significantly from sea state to sea 
state. However, an encounter angle of 90 degrees shows a 
relatively high cost compared to costs calculated for 60 and 
120 degrees at a given sea state. To account for this 
anomaly, the following is suggested. In the regular sea, 
the added mass and inertia were known for a given encounter 
frequency, while in the irregular sea a represent! ve average 
value was used. The method used to obtain the average might 
not represent the actual average. Also, it seems reasonable 
to suppose, that the assumptions of the function weighting 
factor are satisfied for all encounter angles; that is, the 
weighting function (Eg. 3.2), which appears in the cost 
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function (Eg- 3, 1) , does totally represent the added drag 
for all encounter angles. Future study is needed tc answer 
these questions. 

The sea state in the deterministic model is represented 
ty regular waves. On this description, the waves impinging 
on the ship hull ccirespond to only one frequency in the 
energy density spectrum. In the case of irregular seas, 
however, the spectral components change for different 
states, as shown in Figure 6.4. Thus comparison of the 
contrcller parameters obtained for regular seas with results 
for irregular seas is not justified. The function miniriza- 
tion subroutine adjusted contrcller parameters to minimize 
the cost function fcr either case (irregular or regular 
seas) as shown in tables 32 and 33. 



TABLE 29 

SlflOLATIOH RESULTS - STEADY STATE 600 SECS 



FIXED SHIP SPEED (32 KNOTS) IN A REGULAR SEAHAY 
SHIP HGEEL: EQUATIONS OF MOTION 
YAHC=0.0 

ENCOUNTER FREQUENCY = 0.75 RAD/SEC 
ENCOUNTER ANGLE =60 DEGREES 

CONTROLLER A 



sea 

state 

1 

2 

4 

6 



K1 

. 1449795 
.1534657 
. 1514665 
. 1533340 



controller g 
T1 

141.383179 

129.987473 

135.798737 

135.488495 



ains 

T2 

32.9405670 

31.4042358 

32.9749756 

33.5585632 



cost 
J min 

.000764582 

.003056434 

.009345479 

.022174600 



Note that in bcth the deterministic and stochastic 
models, among the controllers considered, contrcller A is 
"BEST” in a seaway disturbance, although the differences in 
cost are not dramatic. 

Finally, the observed dependence of optimal controller 
gains on sea state and encounter angle suggests that an 



48 



TABLE 30 



SIMULATION RESULTS - STEADY STATE 600 SECS 



FIXED SHIP SPEED (32 KNOTS) IN A REGULAR SEAWAY 
SHIP aCDEL: EQUATIONS OF MOTION 
YAHC'O.O 

. ENCOUNTER FREQUENCY = 0-75 EAD/SEC 

SEA STATE 4 



CONTROLLER A 



encounter 




controller qains 


cost 


angle 


K1 


T1 


T2 


J min 


30 


.2351043 


102.021973 


28.3396912 


.00 2985300 


60 


. 1514665 


135.798737 


32.9749756 


.009345479 


90 


. 4964442 


66.546493 


49.7598267 


.048143090 


120 


. 1327230 


149.540543 


33.6013489 


.038937880 


150 


.4536914 


70.566528 


31.5839539 


.062534 153 



TABLE 31 

SIMULATION RESULTS - STEADY STATE 600 SECS 



FIXED SHIP SPEED (32 KNOTS) IN A REGULAR SEAWAY 
SHIP MODEL: EQUATIONS OF MOTION 
YAWC=0.0 

ENCOUNTER FREQUENCY = 0.75 EAD/SEC 
SEA STATE 6 



CONTROLLER A 



encounter 




controller gains 


cost 


angle 


K1 


T1 


T2 


J min 


30 


.2370022 


100.122940 


28.0581207 


.007092878 


60 


. 1533340 


135.488495 


33.5585632 


.022174600 


90 


.5210407 


62.153702 


49,9858093 


. 1 12772880 


120 


. 1414837 


142.695160 


35,3171234 


.091541650 


150 


.4587426 


71.451385 


33.4568024 


, 14461582S 


adaptive controller 


ttust be used 


to provide 


a continue 


minimum on 


the cost 


function. 







After obtaining the optimal gains for controller A, to 
observe the behavior of the rudder and yaw motion of the 
ship, transient response plots were obtained for controller 
A at ship speed of 32 knots and sea state 4 for various 
encounter angles as shown in Figures 6.5 - 6,14. Note the 
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TABLE 32 

CCHPABISCN OF lEBEGOLAfi TO EEGOLAE SEAS COETEOLLEE GAIMS 

SEA STATE 4 



CONTEOILEE A 



encounter 




controller gains 




angle 




K1 


T1 


T2 


30 


(irregular) 


.9815440 


5.733036 


.6999879 


30 


(regular) 


.2351043 


102. 021973 


28.3396912 


60 


(irregular) 


.6201209 


40.805560 


19.6068730 


60 • 


(regular) 


. 1514665 


135. 798737 


32.9749756 


90 


(irregular) 


1.809746 


36.012250 


6.3247080 


90 


(regular) 


,4964442 


66.546493 


49.7598267 


120 


(irregular) 


5. 195190 


18.925130 


. 6999907 


120 


(regular) 


. 1327230 


149.540543 


33.601 3489 


150 


(irregular) 


1.446776 


16. 893750 


. 52654C8C0 


150 


(regular) 


. 4536914 


70.566528 


31. 5839539 



TABLE 33 

COBPAEISCH OF lEEEGDLAE TO EEGOLAE SEAS CONTEOLLEE GAINS 

SEA STATE 6 



CONTBOLLEB A 



encounter 




controller gains 




angle 




K1 


T1 


T2 


30 


(irregular) 


2.9715786 


10-4721832 


. 534 2450 


30 


(regular) 


-23700220 


100. 122940 


28-05812C7 


60 


(irregular) 


1.72-28041 


8-4014740 


.514 1125 


60 


(regular) 


. 15333400 


135.488495 


33.5585632 


90 


(irregular) 


1.8584366 


37. 1672655 


.5792384 


90 


(regular) 


-5210407 


62.153702 


49.9858093 


120 


(irregular) 


3.3422489 


106.722259 


.9260592 


120 


(regular) 


. 1414837 


142.695160 


35.3171234 


150 


(irregular) 


. 2854474 


157. 483887 


119.981018 


150 


(regular) 


. 4587426 


71.451385 


33.4568024 



increase in both rudder and yaw amplitude as the encounter 
angle increased. This is due to the increase in yaw moment 
and sway force on the ship. 
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Pigare 6.5 TAV vs. TIHE 30 DEGREES 
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Figure 6.6 BODDEB vs. TXBE 30 DEGREES 
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Figure 6.7 lAI vs. TIME 60 DEGREES 



\ 
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Figure 6.8 RODDEB vs. TIME 60 DEGREES 
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Figure 6.9 TAB ts. TIHE 90 0E6BSES 



) 



/ 
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Figure 6.10 BDDDEB ys. TIOE 90 DEGBEES 
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32 KNOTS-SEfl STfiTE i 



ENCOUNTER RNGLE 120 DEGREES 




Figure 6.11 YAV vs. TIBE 120 DEGfiEES 



/ 



\ 
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Figure 6.12 EODDEH vs. TIME 120 DEGREES 



r 
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Pigare 6.13 lAl ys* TIHE" 150 DEGBEES 
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Figure 6.14 BODDEE vs. TIME 150 DEGHEES 
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VII- AN ADAP TIVE CONTROLLER 



In a seaway, the controller gains changed dramatically 
for changes in sea state and encounter angle. An adaptive 
controller must be used to provide continuous operation on a 
minimum of the cost function. This Chapter addresses a 
theoretical design of an adaptive controller. 

In the future, there will be better measurement of navi- 
gation than can be provided by conventional equipment on 
board a ship. Presently the Navy is involved in a program 
that will provide precision navigation data. The 
NAVSTAR/GICBAL POSITION SYSTEH (GPS) £Ref. 15] [Ref. 16] 
[fief. 17] will provide extremely accurate three-dimensional 
position and velocity information to users anywhere in the 
world. The position determinations are based on the measure- 
ment of the transit time of EF signals from four satellites 
of a total constellation of eighteen. This system is sched- 
uled to be fully operational in 1988. At present (1984) 
there are four NAVSTAR/GPS satellites in operation which 
allows three to four hours per day of navigation time. 
Already the Texas Instrument Company markets a receiver for 
this system where GPS can be used. 

The Navy Remote Ocean Sensing System (NROSS) [Ref. 18] 
will be able to determine wind velocities over the world's 
oceans with an accuracy sufficient to determine ccean 
surface waves. It's objective will be to acquire global 
ocean data for operation and research use by both the mili- 
tary and civil sectors. This system is scheduled to launch 
its first satellite in June 1989. 

The scheme for an adaptive controller is shown in Figure 
7.1. Having stored the optimal controller parameters in a 
look np table as functions of ship speed, sea state, and 
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encounter angle, the ship operating condition must be Xnown 
so that the table is useful, NAVSTAB/GPS would identify 
ship speed and NBOSS would identify sea state and encounter 
angle. The optimal fararaeters can then be looked up and 
inserted into the controller. This should place system oper- 
ation near the minimun: J. To ensure fine tuning, a micro- 
processor programmed, with the function minimization or-line 
in machine language, with inputs of yaw error and rudder 
motion of the ship would accomplish the fine tuning rapidly. 
Since the subroutine is written in Fortran (as used for this 
study) this would be inappropriate for on-line use. 

The adaptive controller can be performed with digital 
circuits rather than analog components. Garcia [Bef. 19] 
demonstrates the process for converting an analog ccntrciler 
into a digital controller. Figure 7.2 illustrates the 
processing of the major components in a digital controller. 
An analog component circuit can be replaced by an analog to 
digital converter, a digital processor, and a digital to 
analog converter.. Some of the benefits which can be realized 
by doing this are; 

1. A high-speed processor could actually process a 
number of multiplexed signals, performing processing func- 
tions on a number of independent channels. 

2. The processing function is permanent in software, 
unless deliberately changed, and will not drift with age. 

3. The processing function can be changed without 
changing components, merely by changing software. 

4. Accuracy can be made very high and can be changed 
merely by changing software. 

5. Processing, which previously required large compo- 
nents such as inductors in low-frequency controllers, can 
now be performed by very small digital circuits. 
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Figure 7. 1 ADAPTIVE CONTBOILEB 
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Figure 7.2 DIGITAL BLOCK DIAGHAH 

In converting an analog controller to a digital 
controller, the process can he broken down into the 
following steps: 

1. Determine the desired analog transfer function. 

2. Set the sampling frequency. 

3. Apply the bilinear z-transf ormation. 

e- 

4. batch one point in the s domain to the z domain. 

5. Obtain the optimum constant coefficients. 

6. Obtain the digital transfer function. 

7. Obtain the simulation diagram. 

The optimal ccntroller parameters can be stored in 
memory. Intel company markets a 4 megabit non-volatile 
read/ write bubble memory. It is supported by a 7SLI 



t 
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contrcller which provides a black box interface. It is easy 
to use and can be used with any 8- or 16-bit microproces- 
sors. Tfce bubble merocry advantage is: 

1. Fast access time compared with disk or tape. 

2- Ncn-volatile. 

3. Wide temperature range of operation, 

4. Working storage. 

5. Portable operation 

6. Low power, 

?• High reliability. 
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VIII. CO NCLD SIGNS MS BEC OMMENDATION5 FOB FOTORE SIDDY 



A. CONCLUSIONS 

In designing the controller, three different shif acdels 
were used. Using the second-order Nomoto model Equation 2-1 
allowed comparison of results with Reid’s £Bef- 7] £Bef- 10] 
work- It is clear that the answers obtained by function 
minimization agree closely with Reid’s results as shown in 
Tables 4 and 5- A better description of the ship is the 
third-order Nomoto model which involves both the sway and 
yaw equations. This model includes the two dominating poles 
of the ship. The best model to describe the dynamics of the 
ship is a Taylor’s series expansion- This aliowes both 
linear and nonlinear terms in the equations of moticn to 
affect the design of the controller. 

To determine which controller structure would provide 
the minimum cost due to steering, various structures were 
studied. It was found that the dynamics of the plant deter- 
mines the optimum structure for the controller- In calm 
water study, when using a second-order Nomoto model, the 
best structure was controller A. Hhen the third-order Ncmoto 
model Equation 2.2 was used the best structure was 
controller C, but the difference is slight. Observe that in 
each case the controller zeros cancel the plant poles- When 
the equations of motion were used for the plant, the best 
structure was controller C. When the equations cf motion 
were coupled to a sea state generator and the cost function 
was minimized in the presence of a seaway, the best struc- 
ture was ccntroller A. This study concludes that controller 
A should be used. 

A function minimization subroutine is an engineer's tool 
which can be used in many engineering problems. Previously a 
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Batched filter was designed for the Naval Postgraduate 
School research project on the Space Transportation System 
(STS) for the Get Away Special Program. It was matched to 
the signature of the auxiliary power unit (APU) on board the 
space shuttle. The goal was to turn on the solid state 
recording system before lift off, to record the acoustic 
power generated inside the shuttle bay. Basically the 
matched filter is a finite Impulse Response (FIH) filter 
with the weights calculated to obtain the least squared 
error of the desired output when the input is the signature 
of the APO. Figure 8.1 shows the scheme used to evaluate 
the FIR weights. 

B. EICCHHENDATIOMS FCR FOTOEE STUDY 

In the future most ships both military and commercial 
will have GPS receivers as part of their navigation equip- 
ment. Using extremely accurate three-dimensional position 
and velocity information from satellite platforms will allow 
ships to navigate accurately in and out of ports. The func- 
tion minimization subroutine is a powerful tool for 
designing the controller. This routine simply takes the 
inputs that require linimizaticn and adjusts the parameters 
to accomplish this task. The cost function for the added 
drag due to steering is a function of yaw error and rudder 
motion. The use of function minimization and NAVSTAE/GPS 
provides the means for optimization for guidance and cont- 
roll. There are several areas that need future study and 
work . 

1. Should the objective change to track following than 
it is necessary to minimize the yaw error only. This would 
be very important both militarily and commercially should a 
port be mined. If the ship could follow a stringent route, 
knowledge of mine locations would allow access. 
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Finite Impulse Response Filter 




Figure 8.1 HATCHED FILTEB DESIGN 

2. A controller for orbit keeping for satellites with 
High-Energy Laser weapons would be very important. The small 
far- field spot size of a focused laser beam can be selec- 
tively focused on the most vulnerable component on the 
target. Facilitating precision energy deposition, amd 
greatly increasing the probability of a kill. 

3. An adaptive controller to minimize track error on 
board a cruise missile could be programmed for selective 
targets. 
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4. Military and 
do ships ty reducing 



commercial aircraft can benefit just as 
drag to minimize fuel consumption. 
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APPENDIX A 



PROGfiAH 10 CALCULATE OPTIHAL GAINS 

The program is set up to calculate the optimal gains for 
controller A- It is referenced in Chapter five and six. It 
can easily be modified to obtain optimal gains for the rest 
of the ccntrollers. After obtaining the optimal gains the 
program most be modified to do a simulation. The program has 
sufficient comments for appropriate changes. It is refer- 
enced in Chapter two. 

This program can be modified to obtain the Ncmcto 
models. It is referenced in Chapter two. The following need 
to be changed. 

C GAIN COEFFICIENTS TO BE OPTIMIZED 
K=XX{1) 

TP1=XX (2) 

TZ=XX (3) 

TP2=XX (4) 

C EERCB SIGNAL TO EBIVE RUDDER (YAW ACTUAL - YAW CCMMANI) 

C FOR EQUATIONS OF MOTION, 

D=YAW - YAWC 

C ERROR SIGNAL TO DRIVE RUDDER (YAH COMMAND - YAH ACTUAI) 

C FCR NCMOTO 3RD OEIER MODEL. 

E2=YAWC-YAW2 
X1= (D2-X2)/TF1 
X3=K* (TZ + X1 + X2) 

X4= (X3-X5)/TF2 
C INTEGRATION 

X2=X2+X1*DELT 
X5=X5+X4*DELT 
YAW2=YAW24-X5*ESLT 
C COST FUNCTION 

TDIFF=TDIFF+ (YAH-YAW2) **2 
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PB0G3AH TO CALCDIATE OPTIMAL GAINS FOB CONTEOXLEB 



//GAECIA JOB (2220, 0356) , ’RESEARCH' ,CLASS=J 
// EXEC rRTXCLGP,IMSI=DP,REGION=1024K 
//FCRa.SYSIN DD * 

C 

THIS PROGRAM WILL OBTAIN TEE CONTROLLER 
IT IS REFERENCED IN CHAPTER 5- 



OPTIHAL GAINS. 



C 

C 

C 



IN ORDER TO PERFORM 
OBTAINED CHANGE XS(*) 
DIMENSION XS(3 
XS (1)=- 427 
XS(2)=48.66 
XS (3 =10-7 
IS THE STARTING 
IS THE LOWER 
IS THE UPPER 
XL (1) =. 1 
XU 1) = 10. 

XL (2 =1. 

XU 2) =200. 

XL (3) =.10 



SIMULATION ONLY WHEN 



TO X(*) AND DELETE 



C 

C 



c 

c 

c 

c 

c 



) ,XU 



XS (I) 
XL i; 
XU I) 



GUESS 
LIMIT FOR THE 
LIMIT FOR THE 



GAINS 
XU (*) , 



HAVE BEEN 
AND XL (*) 



I'TH 

I'TH 



VARIABLE 

VARIABLE 



A 
IS 



XU (31=100 
DESCfilPTION OF 



THE FOLLOWING 
ECXPLX 



C 

C 

C 



25 

30 

40 



PARAMETERS 

DISCUSSED 
R=9./13- 
NTA=1000 
NPR=100 
NAY=0 
NV=3 
IP=0 

TEE FOLLOWING STATEMENT MUST BE CHANGED TO 
r A T T FT A NT IY\ 

IE ONLY SIMULATICN IS WANTED 

CALL BOXPLX (NV,NAV,NPR,NTA, R,XS,IP, XU,XL,YMN,IER) 
WRITS (6,25) 

■ ' GAINS*,/) 



OPTIMAL 



)=*,F14.7) 



FORMAT hi, * 

DO 30 1=1,3 
KRITE(6,40) I ,XS (I) 

FORMAT (1 X,^X (*,I2; * 

STOP 
END 

SUBROUTINE PLANT (XX) 

SUBHCUTTNE PLANT (XX) SIMULATES THE SHIP 
COMMON TDIFF 
REALMS L,L2,L3,14,L5,L6 

REALMS X, XDOT, 5, YDOT, U , UDOT, V, VDOT , YAW,R,RDOT 
REAL*8 TIME,ETIME,XUDOT,XUU-XVR,XYV,XDD 
BEAL+8 YV, YE,YC,YVVR, YVHR, YVVV, YRRB ,YDDD, YVDOT 
. BEAL*8 NV,NB .NE.NVVR, NVBR, NVVV,NRRE, NDDD, NRDOT 
REAL*8 RHO,I2,FX,FY,MZ,XP,MASS,DELT 
REAX*8 DYAWE, YAWE,YAWC,ISE, ISR,LAMDA, D 
REAL*8 K1,T1, T2,D,X2,DX2,CH ■' " 



DIMENSION XX (3) 



(11J,S 



CLOSE LOOP ANALYSIS WITH FILTER 

INITIAL CONDITIONS FOR INTEGRATION 

SIMULATION END TIME IN SECONDS 
ETIME=600- 
TIME=0- 0 
ICCUNT=1 

INITIALIZE THE COST FUNCTION 
ISE=0.0 
ISR=0-0 
TDIFF=0.0 
LAMDA=4.2 
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C GAIN COEFFICIENTS TO BE OPII.MIZED 
K1=XX{1) 

T1=XX (2) 

T2=XX (3[ 

C X ,XDCI,Y- YDOT ABF FIXED COORDINATES ON EARTH 
X^O-0 

Y = 0.0 
XDCT=0. 0 
YDCI=0. 0 

C D,UDCT, V,VDOT ARE FIXED COORDINATES ON SHIP 

V = 0. 0 
UEOT=0. 0 
VBCT=0. 0 
YAR=0.0 
E=0. 0 
REOT=0. 0 

C CRDEfiEE SPEED IN FEET/SEC 
C 54. Cl FT/SEC=32 KNOTS 
UC=54.01 

C AT STEADY STATE ACTUAL SPEED (U) = COHMAND SPEED (UC) 

U = UC 

C D = RUDDER ANGLE 
D = 0.0 
1=880.5 
12=L**2 
L3=L*L*L 
L4=L*L3 
L5=L*L4 
L6=L*L5 

C SEA DISTURBANCE 

C FCfiCES IN X,Y DIRECTICN COMPUTED IN POUND FORCE 
C HCMENTS IN Z 
FX=0. 

FY=0. 

HZ=0. 

C ISEA IS A SWITCH ; ISEA=0 (CALM WATER) ISEA=1 (SEA STATE) 
ISEA=1 

C HYBRQDYNAMIC COEFFICIENTS ARE INSERTED HERE AS PARAMETERS 
RHC=1.9876 

MASS= (. 0044) * (.5=«RHO^'L3) 

IZ=(0. 00028) * |.5*RHO*L5) 

YAWE=0.0 
X2=0-0 
DX2=0.0 
200 CONTINUE 

S=DSQRT (U**2+V**2) 

C INPUT YAH COMMAND 
YAHC=0.0 

IF (TIME. GE. 0.0)_ YAWC = 0.0 

C ERROR SIGNAL TO EEIVE RUDDER (YAH ACTUAL - YAH ORDERED) 

C ( COMPENSATOR FILTER ) 

YAWE=YAW - YAKC 
DX2= (YAWE-X2) /T2 
D=K1’^ (T1*DX2+X2) 

C AXIAL FORCE HYDRCEYNAMIC COEFFICIENTS (SURGE) 

C XUDOT IS THE ADDFE MASS TERM WHICH MUST BE CHANGED FOB 
C DIFFERENT ENCOUNTER ANGLES , SPEED , ENCOUNTER FREQUENCY 
C 

XUDOT= (-. 000 1) ♦(. 5*RHO*L3) 

XU= (-0. 0253) * (.5*EHO*L2*S) 

XDD= (-0.0 003) * (.5*RHO*L2) 

X7R= (0. 0039) * (.5*RHO*L3) 

XVV= (-. 0012)^ (.5*RHO*L2j 
XDD=(-0.0005) ^ (.5*RHO*L2*S^*2)_ 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 

YV= (-0. 00758) * (. 5*EHO*L2*S) 

YB= (0. 0023)i * (.5*EHO*L3*S) 

YD= (0.001 45) * (.5=»RHO*L2*S**2) 

YVVE=(0.01) *{.5*RHO*L3/S) 
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nn o n on non nnnnnn n nnnnn n nnnnn 



1VEE= 
IVVV= 
YEEE = 
YDDD= 
YODOI IS 
DIPFESEET 



-0. 008) ♦ (.5*RHO=!'L4/S) 
-0.03) ♦ (.5«RHO*12/S) 

0.0031 ’fr .5*RflO*15/S) 

-0. 0005)* (.5*EHO*I2*S=<‘*2) 
THE ADDEE MASS TEEM WHICH 
ENCOUNTER ANGLES, SPEED , 



HOST BE CHANGED FOE 
ENCOUNTER FREQUENCY 



YVDOT= (-0.003 9)* (.5*RHO*L3) 
SPEED=32 KNOTS, ENCCUNTEH 



ANGLE' =150 , ENCOUNTER FREQ =.75 

YVDOT=-2. 3043E+C6 
MCMENT ABOUT Z-AYIS HYDRODYNAMIC COEFFICIENTS (YAH) 

NV= (-0. 00213) *(. 5*RHO *L3*S) 

NB= (-0. 00105) * (,5*EHO*L4*S) 

ND= (-0. 0007) * (.5*RHO*L3*s**2) 

NVVE= (-O.Olo) * (.5*RHO*L4/S) 

NVEE= '-0. 008) *J. 5*RHO*L5/S) 

NVVV= 0.0 1) * (.5*RHO*L3/S) 

-o.oor* ■ • 



NEER= 



)6) * (.5*RHO*L6/S) 



NDDD=J0. 000 1) * (. 5*RHO*L3*S^*2) 

NEDOT IS THE ADDED INERTIA " 

DIFFERENT ENCOUNTER ANGLE 



TERM WHICH MUST BE 
, SPEED , ENCOUNTER 



CHANGED FOR 
FREQUENCY 



NEDOT= (-0.00027) *{.5*RHO*L5) 

SPEED=32 KNOTS. ENCOUNTER ANGLE =150 , ENCOUNTER FREQ =.75 
NRDOT=-1.5096E+1 1 
SETS SEA STATE TO ZERO 

IF ilSEA. EQ. 1) GO TO 30 
FX=0. 

FY=0. 

MZ=0. 

GC TO 35 

TABLE LOOK UP OF SEA DISTURBANCE, 

UNIT 12 HAS THE SEA STATE DATA NAMED CH 
IT MUST BE SYNCHECNIZED BY APPROPRIATELY 
CALLING CH IN THE PROPER TIME IN THE LOOP. 

TEE SEA DATA WAS CREATED FOR 600 SECONDS 
WITH AN INCREMENTAL INTERVAL OF 1 SECOND. 



30 


READ 


(1 




FX= 


CH 


(3 




FY= 


CH 


(4 




MZ= 


CH 


(8 


35 


CCNTII 


JU 



0 ACTUAL SPEED 
UC COMMANDED SPEED 
XF = PROPELLER THEUST 
XP=-XUU*UC**2 
EQUATICNS OF MOTION 

FOR CONSTANT SPEED COMMENT THE NEXT TWO' INSTRUCTIONS 
UDOT= ( (XVR + CASS)*V*E * XUU*U**2 + XVV*V**2 
1 + XDD*D*D + FX + XP ) /(MASS-XUDOT) 

VDOT=(YV*V nR-MASS*D)*R ♦ YD*D + YVVR*V**2*H 

1 + YVBR*V*R**2 + YVVV*V**3 

2 ♦ YERR*R**3 + YDDD*D**3 + FY ) / (MA SS- Y VDOT) 

RDOT= ( NV*V ♦ KR*R ♦ NB*D + NVVR*V**2*R 

1 + NVER*V*R**2 + NVVV*V**3 

2 + NRRR*E**3 + NDDD*D**3 + MZ )/(IZ-NRDOT) 

WEEN TO PRINTOUT 

IF (ICOUNT.EQ.il) GO TO 50 
GC TO 300 

CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW*57.296 
RDEG=E*57.296 
RDDEG=RDOT*57.296 
DDEG=D*57.296 
YAKC=YAWC*57.296 

WRITE (6,100) TIME,XP,X,XDOT, Y, YDOT 
1 , OC,U,UDOT, V,VDOT,YAWC, YAWDEG,RDEG,RDDEG, DDEG 

100 FCRMAT'{1X,*TIME=*,F8. 3, * SEC XP=',F10-2,' LBF X= ' 
1,F8.2,» FT XDCT=» ,F8.4,» FT/SEC Y=',F8.2, 
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n nn noon noonnoooooooonnnn noon onnnno 



2 

3 

4 

5 

6 
7 



II YDOT=* ,F8.4, • FI/SEC* ,/,2x, ' UC=',F8.4, 
FT/SEC 0=’ F8.4,' FT/SEC UD0T= ' ^ FlO . 6 , 

' FT/SEC**2 V=',F8.4,' FT/SEC VDOT=»,F10= 

FT/SEC*^2' ,//2X.’ YAT3C = * - ■ 

DFr: 7AW RATT*=» 



F8.4,' DEG 



7, 



//) 



GO TO 
DEIT 



400 



ON EABTH 



/ t t \J m S f 

i. • ~r 1/AJvj YAW = ' F15« 

DEG YAW BATE='» -F15- 7, ^ DEG^SEC YAW ACCEL=*» 
,115.7,* DEG/£EC**2* ,/,2X, ' EDDDER =',F15.7,» DEG 
ICCUNT=1 

C TEST IF WANT TO STOP 
300 IF (TIME. GE. ETinE) 

C INTEGRATION STEP SIZE 
DELT=1. 0 
C INTEGRATION 

D=0+0DOT*DELT 
V=V+VDOT*DELT 
E=E+RDOT*DELT 
YAW=YAW+R*DEIT 
X2=X2+DX2*DELT 
C CCNVERT SHIP TO FIXED 

XDOT=U*DCOS (YAW) -V*DSIN 
YDOT=U*DSIN (YAW) +V*DCOS 
X=X+XDOT*DELT 
Y=Y+YDOT*DELT 
TIHE=TiaE+DELT 
ICC0NT=IC0UNT+1 
ISE=ISE + LAMIA*YAWE**2 
ISB=ISH + D*»2 
GO TO 200 

J= TDIFF = COST FUNCTION 
TD1FF=ISE+ISE 

WRITS (6-500) ISE,ISE,TDIFF,K 1,T1,T2 

F0EMAT(* ' - 1X, 'ISE=' -F15-7, • ISR=* -F15.7, ' TOTAL=* 
1,F15.7,2x,»K1 = *,F15.7,2X, *T1=* ,F15. 7,2X, *T2 = * ,F15. 7) 
REWIND 12 
RETURN 
END 

DELETE ALL THE FCILOWING SUBROUTINE IF SIMULATICN ONLY 
AND NOT OPTIMIZATION IS WANTED 



COORDINATES 
[YAW) 
YAW) 



400 

500 



SCEEODTINE BOXELX (CATEGORY HO) 

PURPOSE 

BOXPLX IS A SUBROUTINE USED TO SOLVE THE PROBLEM OF 
locacting A MINIMUM (OR MAXIMUM) OF AN ARBITRARY OBJECT- 
iv€ function SUBJECT TO ARBITRARY EXPLICIT AND/OR 
imtlicit constraints by tHE COMPLEX METHOD OF M.J. BOX. 
explicit constraints are dEFINED AS UPPER AND LOWER 
bounds on the independent variables IMPLICIT constraints 
may be arbitrary function of the varlABLES- TWO FUN- 
ction subprogram to evaluate the objective FUNCTION AND 
implicit conSTEAINTS- RESPECTIVELY, must be SUPPLIED 
by the user (see EXAMPLE BELOW). BOXPLX ALSO HAS tHE 
opticn to perform integer programming, where the values 
or the independent variables are restricted to integers. 

USAGE 

CALL BOXPLX (NV , N AV, NPR , NTA, R, XS , IP, XU , XL , YMN ,IER) 



DESCRIPTION OF PARAMETERS 

NV AN INTEGER INPUT DEFINING THE NUMBER OF INDEPENDENT 
VARIABLES OF THE OBJECTIVE FUNCTION TO BE MINIMIZED. 
NOTE: MAXIMUM NV + NAV IS PRESENTLY 50. MAXIMIM NV IS 

25. IF THESE LIMITS MUST BE EXCEEDED, PUNCH A SOURCE 
DECK IN THE USUAL MANNER, AND CHANGE THE DIMENSION 
STATEMENTS. 



74 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



NAV AN INTEGER INEUT DEFINING THE NUMBER OF AUXILIARY var 
iaELES THE USER WISHES TO DEFINE FOR HIS OWN CONVENIENCE- 
TYEICALLY HE MAY WISH TO DEFINE THE VALUE OF EACH IMELICI 
CONSTRAINT FUNCTION AS AN AUXILIARY VARIABLE. IF THIS 
IS DONE, THE OPTICKAL OUTPUT FEATURE OF BOXPLX CAN BE 
USED TO OBSERVE THE VALUES OF THOSE CONSTRAINTS AS THE 
SCIUTICN PROGRESSES. AUXILIARY VARIABLES, IF USED, 

SHOULD BE EVALUATED IN FUNCTION KE (DEFINED BELOW). 

NAV MAY BE ZERO. 



NPB INPUT INTEGER CONTROLLING THE FREQUENCY OF OUTP.UT 
desired for diagncstic purposes. 

IF NPB -LE. 0, NO OUTPUT WILL BE 

PfiOrUCED BY BOXPLX. OTHERWISE, THE CURRENT COMPLEX OF 
K= 2*NV VERTICES AND THEIR CENTROID WILL BE OUTPUT AFTER 
EACH NPR PERMISSIBLE TRIALS. THE NUMBER OF TOTAL TRIALS 
NUMEEE OF FEASIBLE TRIALS, NUMBER OF FUNCTION EVALUAIICN 
AND NUMBER OF IMPLICIT CONSTRAINT EVALUATIONS ARE IN- 
CLUDED IN THE OUTPUT. 

ADDITIONALLY, (WHEN NPR .GT. 0) THE SAME INFORMATION 
HIIL BE OUTPUT: 



1 

I' 

i). 



IF THE INITIAL POINT IS NOT FEASIBLE, 

AFTER THE FIRST COMPLETE COMPLEX IS GENERATED, 

IF A FEASIBLE VERTEX CANNOT BE FOUND AT SOME TRIAL, 
IF THE OBJECTIVE VALUE OF A VERTEX CANNOT BE MADE 
NC-ICNGEE-WORST. 

IF THE LIMIT ON TRIALS (NTA) IS REACHED AND, 

WHEN THE OBJECTIVE FUNCTION HAS BEEN UNCHANGED FOE 
2*NV TRIALS, INDICATING A LOCAL MINIMUM HAS BEEN 
FOUND. 



IF THE USER 
A CHOICE OF 



WISHES TO 
NPR = 25, 



TRACE THE PROGRESS OF A SOLUTION, 
50 OR 100 IS RECOMMENDED. 



NTA INTEGER INPUT OF LIMIT ON THE NUMBER OF TRIALS 
allowed in the calculation. 

IF THE USER INPUTS NTA .LE. 0, A default 
VAIOE OF 2000 IS USED. WHEN THIS LIMIT IS REACHED 
CONTBCL RETURNS TO THE CALLING PROGRAM WITH THE BEST 
ATTAINED OBJECTIVE FUNCTION VALUE IN YMN, AND THE BEST 
ATTAINED SOLUTION POINT IN XS. 

R A REAL NUMBER INPUT TO DEFINE THE FIRST RANDOM NUMBER 
USED IN DEVELOPING THE INITIAL COMPLEX OF 2*NV VERTICIES. 
jo. .GT. R .LT. 1.) IF E IS NOT WITHIN THESE BOUNDS, 

IT WILL EE BEPLACEI ■" ' ^ 



BY 1./3. 



LEAST NV+NAV 



XS INPUT BEAL ARRAY DIMENSIONED AT 
the first nv must contain a 
FEASIBLE ORIGIN FCB STARTING THE CAL- 

COLATICN. THE LASl NAV NEED NOT BE INITIALIZED. UPCN 
RETURN FROM BOXPLX- THE FIRST NV ELEMENTS OF THE ARRAY 
CONTAIN THE COORDINATES OF THE MINIMUM OBJECTIVE 
function, AND THE REMAINING NAV (NAV .GE. 0) CONTAIN THe 
values or THE CORRESPONDING AUXILIARY VARIABLES. 

IP INTEGER INPUT FOR OPTIONAL INTEGER PROGRAMMING, 
if ip=1, THE VALUES OF THE INDEPENDENT VARIABLES WILL 
he replaced WITH INTEGER VALUES (STILL STORED AS REAI*4) . 



XU A REAL 
upper BOUND 
conSTHAiNT) . 

XL A REAL 
lower hound 



ARRAY DIMENSIONED AT LEAST NV INPUTTING THE 
ON EACH INDEPENDENT VARIABLE- (EACH EXPLICIT 
INPUT VALUES ARE SLIGHTLY ALTERED BY BOXPLX 



ARRAY DIMENSIONED AT 
on each independent 



LEAST NV INPUTTING THE 



VARIABLE, (EACH EXPLICIT CONstraint) 
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C NOTE: FOR BOTH XO AND XL CHOOSE REASONABLE 

C VALUES I? NONE ARE GIVEN, NOT VALUES WHICH ARE 
C magnitudes ABOVE CE BELOW THE EXPECTED SOLUTION. 

C input values are SLIGHTLY ALTERED BY BOXPLX. 

C 

C YMN THIS OUTPUT IS THE VALUE (HEAL^U) OF THE OBJECTIVE 
C funcTICN,CO£RESPOKLING TO THE SOLUTION POINT OUTPUT IN XS 
C 

C lEE INTEGER ERROR RETURN. TO BE INTERROGATED UPON 
C return FROM BOXPLX. lER WILL BE ONE OF THE FOLLOWING: 

C 

C =-1 CANNOT FIND FEASIBLE VERTEX OH FEASIBLE CENTROID 
C AT THE START OR A RESTART {SEE 'HETHOD* BELOW). 

C =0 FUNCTION VALUE UNCHANGED FOR 'N* TRIALS. (WHERE 
C N=6*NV+10) THIS IS THE NORMAL RETURN PARAHETER. 

C =1 CANNOT DEVELOP FEASIBLE VERTEX. 

C =2 CANNOT DEVELOP A NO-LCNGEE-WOEST VERTEX. 

C =3 LIMIT ON TRIALS REACHED. (NTA EXCEEDED) 

C NOTE: VALID RESDITS MAY BE RETURNED IN ANY OF THE 

C ABOVE CASES. 

C 

C EXAMPLE OF USAGE 
C 

C THIS EXAMPLE MINIMIZES THE OBJECTIVE FUNCTION SHOWN IN 
C the EXTERNAL FUNCTION FE(X). THERE ARE TWO INDEPENDENT 
C varlAELES X{1) S X (2) , AND TWO IMPLICIT CONSTRAINT 
C function X (3) & X(U) WHICH ARE EVALUATED AS AUXILIARY 
C variatles (see EXTERNAL FUNCTION KE(X) ). 

C 

C DIEENSICN XS (4) ,XD (2) ,XL (2) 

cc 
cc 
c 
c 

cc 
c 
c 

cc 
c 
c 

cc 
c 
c 
c 
c 
c 
c 

cc 
c 
c 
c 

c 

c 

cc 
c 
c 

cc 
cc 
c 

C EVALUATE CONSTRAINTS. SET KE=0 IF NO IMPLICIT CONSTRAINT 
C is viCLATSD, OR SET KE=1 IF ANY IMPLICIT 
c ccnstraint is violated. 

C DIMENSION X (4) 

C XI = X(1) 

C X2 = X{2) 

C KE = 0 

C X(3) = XI + 1.132051 + X2 

C IF (X(3j .LT. 0. .OR. X (3) . GT. 6.) GO TO 1 

C X (4) = X1/1. 732051 -X2 












STARTING GUESS 
XS(1) = 1.0 
XS 2) = 0.5 
UPPER LIMITS 
XU (1) = 6.0 
XU(2) = 6.0 
lOWEH LIMITS 
XI(1) = 0.0 
XLl2) = 0.0 



R = 9./13. 
NTA = 5000 
NPR = 50 
NAV = 2 
NV = 2 
IP = 0 

CALL BOXPLX 



1 

3 ^ 



(NV,NAV,NPE,NTA,R,XS,IP,XU,XL,YMN,IEfi) 
WRITE (6, 1) ( (XS(I) ,1=1,4) ,YMN,IER) 

FORMAT i////, * THE POINT IS LOCATED AT (XS (I) =) ♦ 

, 4 (e 13. / #5x) ,//, 

AND THE FUNCTION VALUE IS *,E13.7,» lER = *,I5) 



STOP 

END 



FUNCTION KE(X) 
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GE. C.) RETURN 



C IF <X(4) . 

CC 

C 1 KE = 1 
C EETUEN 

C END 

CC 
CC 

C FUNCTION FE(X) 

C DIMENSION X(4) 

CC 

CC TEIS IS THE OBJECTIVE FUNCTION. 

C FE= -(X(2)**3 *(9.- (X (1)-3.) **2)/ (46.76538) ) 

C EETUEN 

C END 

C 

C METHOD 

C 

C THE CCMPXEX HETHOE IS AN EXTENSION AND ADAPTION OF 
c tJi€ simple method cf linear programniing . 

C STARTING WITH ANY CNE feasible point in n-dimension 
C A "COMPLEX” Of 2*N vertices is constructed by 
C SELECTING RANDOM FCINTS WITHIN THE feasible 
C HEGICN. FOR THIS PURPOSE N COORDINATES ARE FIRST 
C RANDOMLY CHOSEN WITHIN THE SPACE BOUNDED BY EXPLICIT CCN- 
C STRAINTS. THIS DEFINES A TRIAL INITIAL VERTEX, 
c it is then checked for possible violation 
C OF IMPLICIT CONSTRAINTS. IF one or more are violated, 

C THE TRIAL INITIAL VERTEX IS DISPLACED half of its 
C DISTANCE FROM THE CENTROID OF PREVIOUSLY SELECTED initial 
C VERTICES. IF NECESSARY THIS DISPLACEMENT PROCESS IS 
C REPEATED UNTIL THE VERTEX HAS BECOME FEASIBLE. IF THIS 
c FAIL TO happen after 5*n+10 displacements, 

C THE SCLUTION IS ABANDONED. AFTER EACH VERTEX IS ADDED 
C TC THE COMPLEX, TEE CURRENT centroid is checked fcr 
C FEASIBILITY. IF IT IS INFEASIBLE, the last trail 
C VERTEX IS ABANDONED AND AN EFFORT TO GENERATE an alter- 
C ATIVE TRIAL VERTEX IS MADE. IF 5*N+10 VERTICES ARE 
C ABANDCNED CONSECUTIVELY , THE SOLUTION IS TERMINATED. 

C 

C IF AN INITIAL COMPLEX IS ESTABLISHED, THE BASIC 
c computation loop is initiated. 

C THESE INSTRUCTIONS FIND THE CURRENT WORST vertex, that 
C IS, THE VERTEX WITH THE LARGEST CORRESPONDING value for 
C THE OBJECTIVE FUNCTION, AND REPLACE THAT VERTEX BY 
C ITS OVER-REFLECTICN THROUGH THE CENTROID OF ALL CTHEB 
c vertices, (if the vertex to be 

C REPLACED IS CONSIDERED AS A VECTOR IN n-space, 

C ITS CVER-RSFLECTICK IS OPPOSITE IN DIRECTION, IN- 
C CREASED IN LENGTH BY THE FACTOR 1.3, AND COLLINEAE WITH 
C THE REPLACED VERTEX AND CENTROID OF ALL OTHER VERTICES.) 

C 

C WHEN AN OVER-REFLECTION IS NOT FEASIBLE OR REMAINS 
C WORST, IT IS CONSIBERED NOT-PEEMISSIBLE 
C AND IS DISPLACED HALFWAY TOWARD THE CENTROID. 

C AFTER FOUR SUCH ATTEMPTS ARE MADE UNSUCCESSFULL Y 
C EVERY FIFTH ATTEMPT IS MADE BY REFLECTING THE OFFENDING 
C VERTEX THROUGH THE PRESENT BEST 

C VERTEX, INSTEAD OF THROUGH THE CENTROID. IF 5*n+10 
C DISPLACEMENTS AND OVER-REFLECTIONS OCCUR WITHOUT A 
C SUCCESSFUL (PEHMISSI BL E) RESULT, THE CURRENT BEST 
C VERTEX IS TAKEN AS AN INITIAL FEASIBLE POINT FOR A 
C RESTAET RUN OF THE COMPLETE PROCESS. 

C RESTARTING IS ALSO UNDERTAKEN WHEN 6*nv+10 CONSECUTIVE 
C TRIALS HAVE BEEN MADE WITH NO SIGNIFICANT CHANGE IN THE 
C VALUE OF THE OBJECTIVE FUNCTION. IN ALL CASES, 

C RESTARTING IS INHIBITED IF THE LAST RESTART DID NOT 
C PRODUCE A SIGNIFICANT IMPROVEMENT IN THE MINIMUM 
C ATTAINED. 

C 
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C IT IS EECOMaENDED THAT THE USEE READ THE REFERENCE FOR 
C FURTHER USEFUL INFCEHATION. IT SHOULD BE NOTED THAT THE 
C ALGORITHM DEFINED THERE HAS BEEN ALTERED TO FIND THE 
C CONSTRAINED MINIMUM, RATHER THAN THE MAXIMOM- 
C 
C 

c 

C REMARKS 

C 

C THE INTEGER PROGRAMMING OPTION NAS ADDED TO THIS PROGRAM 
C AS SUGGESTED IN REFERENCE (2) . A MIXED 
c in teqer/continuou£ variable version of boxplx 
C WOULD BE EASY TO CREATE BY DEclarinq *'ip" to be an array 
C OF NV CONTROL VARIABLES WHERE IP (i) = 1 woald indicate 
C THAT THE I-TH VARIABLE IS TO BE CONFINED to integer 
C VALDES* EACH STATEMENT OF THE FORM * IF (IP . EQ. t) ’ etc- 
C WOULD THEN NEED TO BE ALTERED TO »IF (IP (I) - EQ. 1)’ etc. 
C , WHERE THE SUBSCRIPT IS APPROPRIATELY CHOSEN. NORMALLY, 
C XU AND XL VALUES ARE ALTERED TO BE AN EPSILON ’WITHIN' 
c actual values 

C DECLARED BY THE USER. THIS ADJUSTMENT IS NOT MADE 
C WHEN IP=1. 

C 

C NOTE: NO NON-LINEAR PROGRAMMING ALGORITHM CAN GUARANTEE 

c that the answer found is the global 

C MINIMUM, RATHER THAN JUST A local minimum, however, 

C ACCORDING TO REF. 2, THE COMPLEX method has an advantage 
C IN THAT IT TENDS TO FIND THE GLOBAL minimum more 
C FREQUENTLY THAN MANY OTHER NON-LINEAR PROGRAM- 
C MING ALGORITHMS. 

C 

C IT SHOULD BE NOTED THAT THE AUXILIARY YARlkBLZ FEATURE 
c CAN ALSO BE USED TO DEAL WITH 

C PROBLEMS CONTAINING EQUALITY CONstraints. any equality 
C CONSTRAINT IMPLIES THAT A GIVEN VAEiable is not truly 
C INDEPENDENT, THEREFORE, IN GENERAL, ONE variable 
C INVOLVED IN AN EQUALITY CONSTRAINT CAN BE RENUMBERED from 
C THE SET OF NV INDEPENDENT VARIABLES AND ADDED TO THE SET 
C OF NAV AUXILIARY VARIABLES. THIS USUALLY INVOLVES 
C renumbering THE INDEPENDENT VARIABLES OF THE GIVEN 
C problem 

C SUBROUTINES AND FUNCTIONS REQUIRED 
C 

C SUBROUTINE 'BOUT* AND FUNCTION 'FBV* ARE INTEGRAL 
C parts of THE BOXPLX PACKAGE, 

C 

C TWO FUNCTIONS MUST BE SUPPLIED BY THE USER. THE FIRST, 
c ke (x) , is used to evaluate the implicit 

C CONSTRAINTS. SET KE=0 AT THE beginning of the function 
C , TEEN EVALUATE TEE IMPLICIT CONstraints. in the example 
C ABOVE, TEE FIRST CONSTRAINT, X(3j, must be within the 
C RANGE (0. ,LE. X{3) .LE, 6,), THE SECOND constraint x (4) 
C , MUST BE .GE. 0. , IF EITHER CONSTRAINT IS not within 
C THESE BOUNDS, CONTROL IS TRANSFERRED TO STATEMENT 1, 

C AND KE IS SET TO "1" AND CONTROL IS RETURNED TO BOXPIX. 

C 

C THE SECOND FUNCTION THE USER MUST PROVIDE EVALUATES THE 
c objective function, it is 

C CALLED FE{X) AS SHCWN IN THE EXAMple above, and fe 
C MOST EE SET TO THE VALUE OF THE OBJECTIVE function 
C CCBBESPCNDING TO CURRENT VALUES OF THE NV INDEPENDENT 
C VARIABLES IN ARRAY *X'. 

C 

C REFERENCES 

C 

C BOX, M. J., "A NEW METHOD OF CONSTRAINED OPTIMIZATION 
C and a COMPARISON WITH OTHER METHODS", 

C computer journal, 8 apr. '65, PP. 45-52. 
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SOBRCUTINE BOXPLX |NV, NAV , NPR, NTZ ,RZ , XS , IP, BO , BL , YMN ,IER) 

DIMENSION V (50,50), FUN (50) , SUM (25), CEN (25) , XS(NV) 
1 ,hu (nv) , bl (NV) 

KV = 5 
EF = 1.E-5 
NTA = 2000 

IF (NTZ.GT.O) MA = NTZ 
~ HZ 

IF (R.LE.O. .OR.R.GE. 1.) R=1./3. 

NVT = NV+NAV 

TOTAL VARS, EXPLICIT PLUS IMPLICIT 

NT = 0 

CURRENT TRIAL HO. 

NET = 0 

CURRENT NO. OF PERMISSIBLE TRIALS 
NTFS = 0 

CURRENT NO. OF TIMES F HAS BEEN ALMOST UNCHANGEE 
CHECK FEASIBILITY OF START POINT 



2 

•3 



DO 

VT 

IF 

II 

VI 

GC 

IF 

II 

VT 

IF 



4 1=1, NV 
= XS(I) 
JBL^I) .LE 

= BL(I) 

TO 2 

(BO (I) .GE 



= BU(I)i 
(NPR.GT. 



VT) GO TO 1 



VT) GO TO 3 



VJI,1) = VT 
CENa) = VT 
IF ilP. EQ. 

BI (I) = B1(I) ♦ 



0) KRITE (6,49) II 



4 SOM (I) 



GC TO 4 

♦ AMAX1 (EP,EP=«‘ABS (BL (I) ) ) 
AMAX1 (EP,EP*ABS (BU (I) ) ) 



NCE = 1 

NUMBER CF CONSTRAINT EVALUATIONS 

I = 1 
IF 
IF 

WRITE 
GC TO 
NFE = 1 



(KE(V (1, 1) ) .EQ. 0) GO TO 5 
JNPR.LE.O) go TO 12 
- ^ 6 , 50 ) 



RUKEEH OF VERTICES (K) = 2 TIMES NO. OF VARIABLES. 
K = 2*NV 



NUMBER OF DISPL ACEMENTS ALLOWED. 
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c 

c 

c 



c 

c 



c 

c 



c 



c 

c 

c 

c 

c 



c 

c 



c 

c 



c 

c 

c 



KLIi'l = 5^NV+10 

NUaEEB OF CONSECUTIVE TRIALS WITH UNCHANGED FE TC 
ter nirate. 

NCT = NLIH+NV 
AIPHA = 1-3 
FK = K 
FKM = FK-1- 
BITA = ALPHA+1. 

INSURE SEED OF RANDOM NUMBER GENERATOR IS ODD. 

ICR = R*1.E7 

IF (MOD (IQR,2) . EQ.O) ICR=IQR+101 

SET DP INITIAL VERTICES 
FUN{1) = FE(V(1,1)) 

YKN = FUN (1) 

6 FI = 1. 

FUNCLD = FDN ( 1) 

DC 15 1=2, K 
FI = FI+1. 

LIHT = 0 

7 LIMT = LIMT+1 

END CALCULATION IF FEASIBLE CENTROID CANNOT BE FOUND. 
IF (IIMT. GE. Nllfl) GO TO 11 

DC 8 J=1,NV 

RANDOM NUMBER GEKEEATOE (EANDU) 

IQR = IQR*65539 

IF (IQR.LT.O) IQR = IQR+2 147483647+ 1 
ECX = IQR 

EQX = RQX*. 4656613E-9 

VJJ,I) = BL (J) +RQX* (BU ( J)-BL (J) ) 

IF (IP-EQ.1) V (J,I)=AINT (V {J,I) +- 5) 

8 CONTINUE 

DO 10 L=1,NLIM 
NCE = NCE+1 

IF (KE(V{1,I) ).EQ.O) GO TO 13 



DO 9 J=1,NV 
VI = .5* (V (J 
IF (IP.EQ.1) 
V (J,I) = VT 
S CCNTINUE 

10 CCNTINUE 



I)+CEN(J1 ) 

= AINT (VT + .5) 



11 IF (NPR-LE.O) GO TO 12 
HBITE (6,51) I 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,I,FUN,CEN,I) 

12 lER = -1 
GC TO 48 

13 DO 14 J=1,NV 

SUM (J) = SUM (J) +V (J,I) 

14 CEN (J) = SUM ( J)/FI 

TRY TO ASSURE FEASIBLE CENTROID FOR STARTING. 

NCE = NCE+1 

IF (KE (CEN) .EQ.O) GO TO 60 
SUMJj) = SUM(J) -V(J,I) 

GC TO 7 

6C NEE = NFE+1 

FUN (I) = FE(V(1,I)) 

15 CCNTINUE ^ ^ 



I 
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END CF LOOP SETTING OF INITIAL CO^IFLEX. 

IF (NPR.LE.O) GO TO 17 

CALL BOUT (NT,NPT,NFE,NCE, NV, NVT,V, K,FUN,CEN,0) 

FIND THE WORST VERTEX, TEE *J»TH. 

J = 1 

DC 16 1=2, K 

IF_ (FUN (J) . GE.ION (I) ) GGTO 16 

16 CCNTINOE 

BASIC LOOP. ELIMINATE EACH WORST VERTEX IN TURN, 
it Bust Lecome NC LONGER HORST, NOT MERELY IMPROVED, 
find next-to-vertex, THE *JN*TH ONE. 

17 JN = 1 

IF (J.EQ. 1) JN = 2 

DC 18 1=1, K 

IF (I.EQ.J) GC TO 18 

IF (FUN (JN) .GE.FUN (I) ) GO TO 18 

JN = I 

18 CONTINUE 

II«T = NUMBER OF MOVES CUEING THIS TRIAL TOWARD THE 
centxoid DUE TO FUNCTION VALUE. 

LIMT = 1 

CCKPUTE CENTROID AND OVER REFLECT WORST VERTEX. 

DO 19 1=1 ,NV 
VI = V(I,J) 

SUM (I) = SUM(I)-VT 

CEN (ij = SDM(lj/FKM 

VT = BETA+CEN (I) -ALPHA*VT 

IF (IP.EQ.1) \I = AINT(VT + .5) 

INSURE THE EXPLICIT CONSTRAINTS ARE OBSERVED, 

19 V(I,J) = AMAX1 (AMINl (VT,BU (I) ) ,BL (I) ) 

NT = NT+1 

CHECK FOR IMPLICIT CONSTRAINT VIOLATION. 

20 DC 25 N=1,NLIH 
NCE = NCE+1 

IF (KE (V (1, J) ) .EQ, 0) GO TO 26 

EVERY *KV*TH TIME, OVEE-EEFLECT THE OFFENDING VERTEX 
through the BEST VERTEX. 

IF t«OI> (N,KV) . KE.O) GO TO 22 
CALL FBV (K,FUN,M) 



DC 21 1=1, NV 

VT = BETA*V (I,C) -ALPHA + V (I, J) 

IF (IP.EQ.1) \1 = AINTiVT+.5) 

1 V(I,J) = AMAX 1 (AMINl (VI,BU (I) ) , 



BL(I)) 



GC TO 24 

CCNSTEAINT VIOLATION: 



MOVE NEW POINT TOWARD CENTROID. 



22 DC 23 1=1 ,NV 

VT = .5* (CEN (I) +V (I, J) ) 

IF (IP.EQ.1) 1(1= AINT(VT + .5) 



23 



(IP.EQ. 1) 
V(I,J) = VT 
CCNTINUE 
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24 NT = NT-H 

25 CONTINUE 



c 

C 

C 



C 

c 

c 

c 



c 

c 



c 

c 

c 

c 



HR = 1 

CANNOT GET FEASIBLE VERTEX BY MOVING TOHARD CENTROID, 
OR BY OVER-REFLiCTING THRU THE BEST VERTEX. 

IF (NPR.LE.O) GO TO 42 
WRITE (6,52) NT,J 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V, K,FUN,CEN, J) 

GC TO 42 

FEASIBLE VERTEX FOUND, EVALUATE THE OBJECTIVE FUNCTICN. 

26 NFE = NFE+1 

FUNTRY = FE (V (1, J) ) ^ 

TEST TO SEE IF FUNCTION VALUE HAS NOT CHANGED, 

AFO = ABS (FONTEY-FUNOLD) 

AMX = AW AX1 (AES (EP*FUNOID) ,EP) 

ACTIVATE THE FOLLOWING TWG STATEMENTS FOR DIAGNOSTIC 
purposes only. 

WRITE (6,99) J,AFO,AMX,FUNTRY,FUNOLD,FUN (J) ,FUN (JN) 

1 ,NTFS,N 

99 FORMAT ( 1 X, 13 ,6E 1 5. 7 , 215) 

IF (AFO. GT. AMX) GO TO 27 
NTFS = NTFS+1 
IF (NTFS.LT. NCT) GO TO 28 
lER = 0 

IF (NPR.LE.O) GO TO 42 
WRITE (6,53) K 

CALL BOUT (NT,NPT,NFE,NCE,NV, NVT,V,K,FUN,CEN,0) 

GO TO 42 

27 NTFS = 0 

IS TBE NEW VERTEX NO LCNGEE WORST? 

28 IF (FUNTBY.lt. FUN (JN) ) GO TO 34 

TRIAL VERTEX IS STILL WORST: ADJUST TOWARD CENTROID. 
EVERY *KV*TH TIFF. OVER-REFLeCT THE OFFENDING VERTEX 
tircuqh the BEST VERTEX. 

LIWT = LIMT+1 

IF (HOD (LIMT, KV) .NE. 0) GO TO 30 
CALL FBV (K,FDF,M) 

DC 29 1=1 ,NV 

VT = BETA*V (I ,H) -ALPHA*V (I,J) 

IF (IP, EQ. 1) VI = AINT(VT + .5 

29 V(I,J) = AMAX 1 (AMIN1 (VT,BU (I) ) ,BL (I) ) 

GC TO 32 



30 DO 31 1=1. NV 

VT = .5* (CEN (I)+V (I, J) ) 
IF (IP.EQ. 1) VT = AINT I 
VJI.J) = VT 

31 CCNIINUE 



(VT+.5) 



32 IF (LIMT.LT.NIIM) GO TO 33 

CANNOT MAKE THE »J*TH VERTEX NO LONGER WORST BY 
displacing toward 

THE CENTROID OR BY OVER-REFLECTING THRU THE BEST VERTEX. 
lER = 2 

IF (NPR .LE. 0) GO TO 42 
{6 52) ftT J 

CALL BOU’t (NT,NPi,NFE,NCE, NV, NVT,V, K,FUN,CEN, J) 

GC TO 42 

33 NT = NT+1 
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GC TO 20 



SUCCESS: WE HAVE A EEPLACEaENT FOE VERTEX J. 

34 FON(J) = FUNTBY 
FUNOLD = FUNTEY 
NPT = NPTH“1 

EVERY 100»TH PERMISSIBLE TRIAL, RECOMPUTE CENTROID 
suffimation to AVCID CREEPING ERROR. 

IF (MOD (NPT, ICC) -NE- 0) GO TO 37 

DO 36 1=1, NV 
SOM (I) = 0. 

DC 35 N=1 ,K 

35 SUM (I) = SOM (I) +V (I,N) 

CEN(I) = S0H(I)/FK 

36 CCNTINUE 

LC = 0 
GO TO 39 

37 DO 38 1=1, NV 

38 SOM (I) = SUM (I)+V (I, J) 

LC = J 

39 IF (NPR.LE.O) GO TO 40 

IF (MOD (NPT,NPB) .NE.O) GO TO 40 

CALL BOUT {NT,N?T,NFE,NCS,NV,NVT,V, K,FUN,CEN,IC) 

EAS THE MAX. NOMEER OF TRIALS BEEN REACHED HITHCUT 
convergence? IF KOT, GO TO NEW TRIAL. 

40 IF (NT.GE.NTA) GO TO 41 

NEXT-TC-WORST VERTEX NOW BECOMES WORST. 

J = JN 
GC TO 17 

41 lER = 3 

IF (NPR.GT.O) WRITE (6,54) 



COIIECTOR POINT FOE ALL ENDINGS. 



1 
2 
3 

§1 



CANNOT DEVELCE FEASIBLE VERTEX. 

CANNOT DEVELOP A NO- LCNGER-WORST VERTEX. 
FUNCTION VALUE UNCHANGED FOR K TRIALS. 
LIMIT ON TRIALS REACHED. 

CANNOT FIND FEASIBLE VERTEX AT START. 



lEE = 1 
lER = 2 
lEE = 0 
lER = 3 
lER = -1 



42 CCNTINUE 

FIND BEST VERTEX. 

CALL FBV (K,FUF,M) 

IF (IER.GE.3) GO TO 44 

RESTART IF THIS SOLUTION IS SIGNIFICA NTLY BETTER THAN 
the previous, OR IF THIS IS THE FIRST TRY, 

IF (NPR.LE.O) GO TO 43 
WRITE (6,55) (P-YMN, FUN (M) I 

43 IF (FUN (M) .GE.YMN) GO TO 47 

IF (ABS {FUN(M)-YMN).LE.AMAX1 (EP,EP*YMN)) GO TO 47 

GIVE IT ANOTHER TRY UNLESS LIMIT ON TRIALS REACHED, 

44 YMN = FUN(M) 

FUN(1) = FON(M) 

DO 45 1=1, NV 
CIN (I) = V(I,M) 

SDM(I) = V(I/M) 
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C 



C 



C 



C 



C 



c 



c 

c 



c 



45 Vjl,1) = V(I,K) 

DC 46 1=1, NVT 

46 XS(I) = V(I,M) 

If (IER.LT.3) GO TO 6 

47 If (NPR.LE.Ol GO TO 48 

CAll BOUT (NT,NPT,NFE,iNCE,NV,NVT,V,K,fUxN,V (1,M) ,-1) 
WRITE (6,55) fUN(M) 

48 RETURN 

49 FORMAT (50H0INDEX AND DIRECTION OF OUTLYING 
1 variable at startiS) 

50 FORMAT (SO HO IMPLICIT CONSTRAINT VIOLATED AT 
Istart. dead eud.l 

51 FCEMAT ('OCANNCT FIND FEASIBLE *, 14, *TH VERTEX OR 
Icentroid at start.*) 

52 FCHHAT (10H0AI TRIAL I4,54H CANNOT FIND FEASIBLE 
Ivertex which is no 

2LCNGER WORST, 14, 15X, 'RESTART FROM BEST VERTEX.*) 

3 FORMAT i40H0FUNCTION HAS BEEN ALMOST UNCHANGED 
Ifcr i5,7h trails) 

4 FORMAT (27H0LIf'.IT ON TRIALS EXCEEDED. ) 

5 FCRMAT(*03EST VERTEX IS NO.*,I3,*OLD MIN WAS*, £15. 7, 

1 * NEw MIN IS * E15. 7) 

56 FORMAT (»0MIN CB^ECTIVE FUNCTION IS »,E15.7) 

END 

SUBROUTINE FBV (K,FUN,M) 

DIMENSION FUN (50) 

2=1 



DC 1 1=2, K 
IF_ (FUN (M) , 

CONTINUE 

RETURN 

END 

SUBROUTINE 



LE.IUN (I) ) GO TO 1 



BOOT 1NT,NPT,NFE, NCE.NV, NVT,y,K,FN,C, IK) 
DIMENSION V(50,50), FN(50), C (25) 

WRITE (6,4) NT,NPT,NFE,NCE 



DC 1 1=1, K 
WRITE (6,5) 
IF (NVT.LE. 
NVP = NV+1 
WRITE (6,6) 
CONTINUE 



FN 

NV) 



(V(J,I) ,J=1,NV) 
GD TO 1 



(V (J,I) , J=NVE,NVT) 



IF (IK.NE.O) GO TO 2 

WRITS (6,7) (C(I) ,I=1,NV) 

RETURN 

2 IF (IK.GE.O) GO TO 3 

WEIJE (6,8) (C(I),I=1,NV) 

E E T D R N 

3 WRITE (6,9) IK,(C(I) ,I=1,NV) 

RETURN 

4 FORMAT (*ONO. TOTAL TRIALS = *,I5,4X, 

1*no- feasible trails = *,i5,4x, 

2*NO. FUNCTION EVALUATIONS = *,I5,4X, 

3*nc. constraint evaluations = *,i5/ 

4*0 FUNCTION VALUE * ,6X ,* INDEPENDENT VARIABLES/ 

SdependENT OR IMPLICIT CONSTRAINTS *) 

5 format (1H ,E18.7,2X,7E14.7/(21X,7E 

6 FORMAT '21Xe7E14.7) 

7 FCEMAT (lOHOCENTEOxD 1 1 X, 7 El 4 .7/ (2 1 X . 7E1 4. 7) ) 

8 FORMAT (*0 BEST VERT EX * , 7X, 7E1 4 . 7/ (2 IX, 7 E 1 4 . 7) ) 



14.7)) 
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9 FORMAT (*OCEN1EOID LESS V X ' , 12 , 2X , 7 E 1 4 
17e14.7) ) 

END 

FUNCTION FE(X) 

DIMENSION Xj3) 

CCMMCN TDIFF 
CALL PLANT (X) 

FE=TDIFF 

RETURN 

END 

FUNCTION KE(X) 

DIMENSION X(3) 

KE=0 

RETURN 

END 

//GO.SYSIN DD * 

//GO. FT12F001 DD DISP=SHR,DSN=MSS. S2160. A34 1 



7/(21X, 
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A20NDIX B 

EXAflPIE PROBLEM OSIHG ICSOS 



The purpose of this example is to demonstrate the 
performance of the program. Consider the control system of 
Figure 4.1 with controller C. Figure B. 1 shows the fclock 
diagram to evaluate the controller parameters. 
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Figure B. 1 DETAIL BLOCK OIAGRAH 





Ih€ differential equations describing the system and its 
desired performance are: 

22=DX2 
Xi4=CX4 
X5=DX5 
E =DR 

ylw=E 

Defining the following cost function: 

J= y'(LAI1DA*YAWi**2 + D**2) dt 
Defining the special functions: 

YAHE=YAWC-YAW 
DJ2= (YAWE-X2) /12 
X3=K1* (T1*DX2+X2) 

DX4= (X3-X4) /T4 
d= (t3"*dx4+x4) 
dx5= (d-x5)/tp1 
x8=k* (tz*dx5 + x5) 
dr= (x8-r) /tp2 
Defining the constants: 

YAHC=1. 0 
K=. 14771 
12=11.2833 
TF1=6.4699 
TF2=53,7931 
1AMDA=4. 2 

and using YAWC=1-0 the optimal solution found by 
the program is: 

K1=. 4179916 
11=53.69932 
12=4.970023 
13=6-294369 
14=13.85919 
CCS1=68. 04735 



88 



Table 34 shows the specifications of this problem with 
the free parameter optimum values found. Figure B.2 shows 
the actual yaw and rudder response. 
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TABLE 34 
ICSOS OOTPDT 



SPECIFICATICNS 

VARIABLES £ INITIAL CCNCITIONS^ 

X2 = .0 

X4 = .0 

X5 = .0 

R = ,0 

YAW = •€ 

J = oO 
TIME = *0 

CONSTANTS: 

YAWC = loOOOOOCOOO 
K = •1477100000 
TZ = 11.283300CC 
TPl = 6.469900COO 
TP2 = 52o793100C0 
LAMDA = 4.20COCCOOO 



FREE PARAMETERS: 



Kl 


CV = 


.4179916 


LL= clCOCOOO 


UL = 


1.000030 


Tl 


0V= 


53.69928 


LL- l.COOOOO 


LL = 


ICO. 000 0 


T2 


OV- 


4.970C29 


LL= l.OOCOOO 


UL= 


!C. 0030 0 


T3 


ov= 


6.294369 


LL= i.oocnoo 


UL = 


10.00000 


T4 


ov= 


13. 85917 


LL= l.COGOOO 


LL = 


2C.0000C 



SPECIAL FUNCTICNS: 

YAWE = YAWC- YAW 
DX2 = (YAWE-X2)/T2 
X3 = K1«(T1*0X24X2 ) 

DX4 = (X2-X4)/T4 
D * (T3>CX4+X4) 

DX5 = (C-X5)/TFl 
' X8 = K*(TZ*D X54)(5) 

DR = (X8-R)/TP2 

DERIVATIVES: 

D(X2 /D(TIME ) = = 

DX2 

D(X4 /D(TIM£ ) = = 

DX4 ^ 

DJX5 /D(TIME ) = = 

DX5 

D(R /D{TIME ) = = 

OR 

D(YAW /D(TIME ) = = 

R 

D(J /D(TIME ) = = 

LAKDA*YAWE**2+D’>=»2 



OUTPUTS: 

TITLE: ACTUAL YAW 
TABULATE: TIFF. 

AT INTERVAL 
PLOT: 0 YAW 

AGAINST: TIME 

END' CALCULATION WHEN 



V 



AND RUDDER RESFCNSE 
D R YAW 

2oCC0C 00 0 0 0 

AT INTERVAL 2.000 OOOCOO 
TIME oGEo 600. OCO 
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Figure B,2 YAH AND RODDEB vs. TIME 
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